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Abstract 

This  thesis  is  the  latest  extension  of  a  line  of  re¬ 
search  begun  over  ten  years  ago.  The  purpose  of  this  re¬ 
search  has  been  to  develop  an  algorithm  to  track  airborne 
targets  (aircraft  and/or  ballistic  missiles)  using  forward 
looking  infrared  (FLIR)  measurements,  as  a  means  of  aiming  a 
high  energy  laser. 

This  research  deviates  from  past  research  in  considering 
the  use  of  actively  acquired  measurements.  Past  research 
has  concentrated  on  the  use  of  passively  acquired  measure¬ 
ments,  i.e.  measurements  of  the  target's  thermal  intensity 
functions  (hotspots)  from  an  array  of  infrared  detector 
elements.  This  research  considered  illumination  of  the 
target  by  a  low  power  laser.  The  measurement  of  the  re¬ 
flected  laser  light  would  then  give  information  about  the 
hardbody  location,  and  presumably  an  aiming  point  for  the 
high  power  laser. 

Specifically,  this  thesis  investigated  a  ballistic 
missile  in  boost  phase  of  flight.  Measurements  of  the 
missile  exhaust  plume  thermal  intensity  from  an  array  of 
infrared  detector  elements  were  used  by  an  enhanced  cor¬ 
relator/linear  Kalman  filter  to  produce  estimates  of  the 
FLIR  image  centroid  location  and  velocity.  These  estimates 
were  then  used  to  simulate  the  aiming  of  a  low  power  laser 
at  the  missile  hardbody.  The  "pseudo-measurement"  output  of 
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an  optical  sensor  receiving  the  reflections  from  the  missile 
hardbody  was  then  used  by  a  second  Kalman  filter  to  estimate 
the  location  of  the  missile  mass  center. 

This  thesis  effort  involved  sensitivity  and  robustness 
studies  of  the  measurement  noise  variance  in  the  filter 
which  estimates  the  missile  mass  center.  These  studies 
indicated  the  filter's  relative  insensitivity  to  changes  in 
the  measurement  noise  variance;  this  parameter  only  affected 
the  transient  time  for  the  filter  to  reach  the  steady  state 
value  of  the  mass  center  location.  Other  parameter  studies 
were  conducted  involving  the  distance  between  the  missile 
mass  center  and  the  exhaust  plume  intensity  center,  and  the 
infrared  sensor  element  resolution  size.  The  first  study 
indicated  decreased  filter  performance  in  locating  the  mass 
center  with  increased  distance.  The  results  from  the  second 
study  were  inconclusive  and  require  further  work. 
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ENHANCED  TRACKING  OF  BALLISTIC  TARGETS  USING 


FORWARD  LOOKING  INFRARED  MEASUREMENTS 
WITH  ACTIVE  TARGET  ILLUMINATION 


1 .  Introduction 


With  the  advent  of  the  Strategic  Defense  Initiative 
(SDI),  research  into  the  use  of  the  high  energy  laser  as  a 
weapon  has  intensified.  The  laser's  ability  to  concentrate 
energy  in  a  small  area  makes  it  especially  attractive  for 
use  against  airborne  targets  (aircraft  and/or  ballistic  mis¬ 
siles).  Any  system  using  a  high  energy  laser  must  have  the 
capabilities  of  autonomously  tracking  targets  and  accurately 
pointing  the  laser  at  those  targets.  It  is  the  requirement 
for  autonomous  tracking  that  motivates  this  line  of 
research . 

1 . 1  Background 

The  Air  Force  Weapons  Laboratory  ( AFWL )  at  Kirtland  Air 
Force  Base,  New  Mexico,  is  presently  engaged  in  research 
involving  the  use  of  high  energy  lasers  against  airborne 
targets.  In  the  system  now  under  development,  targets  are 
passively  detected  using  a  forward  looking  infrared  ( FLIR ) 
sensor,  consisting  of  a  300  x  500  pixel  array  of  detectors. 
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Each  pixel,  or  picture  element,  detects  infrared  energy  over 
an  area  15  microradians  (azimuth)  by  15  microradians 
(elevation)  square  f 1 4  ]  .  For  this  and  previous  work 
[13,14],  a  subset  of  this  array  consisting  of  an  8  x  8  array 
of  pixels  is  used  to  process  infrared  energy  for  target 
tracking.  This  subarray  of  pixels  is  defined  as  the  track¬ 
ing  field  of  view  (FOV)  [13,14]. 

In  the  system  under  consideration,  the  laser  and  FLIR 
sensor  share  a  common  optical  aperture.  The  tracking 
algorithm  uses  the  incoming  FLIR  data  to  determine  a 
position  offset  between  the  target  position  and  the  current 
center  of  the  FOV.  A  controller  then  zeros  these  offsets  to 
keep  the  center  of  the  sensor  FOV  on  the  target.  Because  of 
the  shared  aperture,  centering  the  sensor  FOV  on  the  target 
also  insures  the  laser  is  pointing  at  the  target. 

In  earlier  research,  the  tracking  function  was  accom¬ 
plished  via  a  standard  correlation  algorithm.  This  algo¬ 
rithm  compares  the  current  target  FLIR  image  with  target 
image  data  from  the  previous  sample  period.  The  correlation 
between  past  and  present  image  data  is  used  to  generate  the 
relative  target  position  offsets,  since  a  translation  of  the 
target  image  in  the  FLIR  sensor  plane  is  assumed  to  repre¬ 
sent  a  spatial  translation  of  the  actual  target.  This 
algorithm  has  the  advantages  of  relative  real-time  implemen¬ 
tation  ease,  and  of  having  relatively  good  tracking  perfor¬ 
mance  against  a  wide  range  of  targets,  particularly 
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spatially  distributed  ones.  However,  the  use  of  a 
correlation  algorithm  does  have  some  associated  disad¬ 
vantages  . 

First,  an  inherent  time  lag  exists.  This  lag  is  due  to 
the  time  required  for  computation  of  target  image  correl¬ 
ation  as  well  as  the  lag  associated  with  physical  pointing 
of  the  tracker  at  the  target.  Second,  the  correlation  algo¬ 
rithm  is  unable  to  distinguish  between  image  spatial  trans¬ 
lation  due  to  target  dynamics  and  image  translation  due  uo 
other  factors  such  as  atmospheric  jitter  [10]  and  platform 
bending/vibration  [4].  Finally,  the  correlation  algorithm 
is  unable  to  take  advantage  of  information  on  target  charac¬ 
teristics,  such  as  size  or  shape,  which  may  be  known  a 
priori  or  could  be  estimated  in  real  time. 

To  overcome  these  disadvantages,  Kalman  filtering 
methodology  has  been  incorporated  into  the  tracking  algo¬ 
rithm  for  the  purposes  of  estimation  and  control  [1-5,11- 
18].  By  modeling  the  effects  of  target  dynamics,  atmo¬ 
spheric  jitter,  and  platform  bending/vibration,  and  includ¬ 
ing  them  in  the  filter  dynamics  model,  enhanced  estimates  of 
target  position  can  be  produced.  These  estimates  can  then 
be  propagated  forward  to  produce  a  position  estimate  at  some 
future  time  for  use  in  laser  pointing  and  target  tracking, 
thereby  providing  performance  potential  superior  to  that 
accomplishable  by  a  simple  correlation  algorithm. 
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1 . 2  Summary  of  Previous  AFIT  Research 

Since  1978,  The  Air  Force  Institute  Of  Technology 
(AFIT)  has  had  a  continuing  line  of  research  investigating 
the  use  of  Kalman  filtering  techniques  in  the  AFWL  high 
energy  laser  pointing  and  tracking  system.  This  research 
has  produced  numerous  theses  and  papers.  Previous  theses 
[4,13,14]  have  summarized  that  body  of  work.  That  summary 
is  included  here  with  modifications. 

The  first  study,  begun  in  1978  by  Mercier  [10],  com¬ 
pared  an  extended  Kalman  filter  (EKF)  with  a  standard  cor¬ 
relator  tracker.  The  targets  were  represented  as  point 
sources  of  infrared  energy,  and  the  FLIR  sensor  image  was 
modeled  as  a  Gaussian  bivariate  distribution  with  circular 
equal-intensity  contours.  The  four-state  filter  maintained 
estimates  of  target  position  and  atmospheric  jitter  states 
in  each  of  two  FLIR  plane  coordinate  directions.  Both 
position  and  jitter  states  were  modeled  as  first  order, 
zero-mean,  Gauss-Markov  processes.  For  that  effort,  FLIR 
measurement  noise  was  modeled  as  consisting  of  background 
clutter  and  internal  sensor  noise  (thermal  noise  and  dark 
current).  The  measurement  noise  was  also  considered  to  be 
temporally  and  spatially  uncorrelated.  The  results  showed 
the  EKF's  ability  to  outperform  the  standard  correlator  in 
tracking  long-range  benign  trajectory  targets  by  about  an 
order  of  magnitude  in  rms  tracking  error. 
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Harnly  and  Jensen  [2]  expanded  Mercier's  work  by  ex¬ 
amining  target  scenarios  under  less  restrictive  assumptions 
Target  velocity  and  acceleration  states  were  added  to  the 
filter  in  order  to  handle  less  benign  target  trajectories, 
and  the  FLIR  plane  target  image  was  modeled  with  elliptical 
rather  than  circular  egual-intensity  contours.  Adaptive 
estimation  of  parameters  associated  with  the  target  shape 
function  was  utilized  to  provide  desired  performance.  The 
measurement  noise  model  was  also  modified  to  allow  for 
spatially  correlated  background  noise.  Adaptive  estimation 
of  filter  dynamics  driving  noise  was  also  incorporated  into 
the  tracking  algorithm.  The  results  indicated  the  filter's 
increased  ability  to  track  targets  over  a  wider  range  of 
maneuver  scenarios. 

The  research  conducted  by  Mercier,  and  Harnly  and 
Jensen  had  been  based  on  the  assumption  that  the  target 
image  intensity  function  was  known  a  priori  to  be  composed 
of  a  single  hot  spot  and  well  modeled  as  a  bivariate  Gaus¬ 
sian  function.  Singletery  [16]  and  Rogers  [15]  pursued  a 
line  of  research  which  made  no  such  assumption.  The  algo¬ 
rithms  developed  used  exponential  smoothing  of  centered 
target  images,  as  an  approximation  to  finite  memory  averag¬ 
ing,  to  produce  an  estimate  of  the  target  image  intensity 
shape  function,  or  template.  This  shape  function  was  then 
used  directly  in  the  measurement  model  of  an  Kalman  filter. 
These  algorithms  were  then  tested  using  target  scenarios  in 
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which  targets  were  represented  as  multiple  hot  spots,  and 
the  results  indicated  performance  comparable  to  that  of 
filters  used  previously  on  single  hot-spot  targets,  in  which 
only  parameters  in  an  assumed  bivariate  Gaussian  intensity 
had  to  be  identified  rather  than  an  entire  intensity  func¬ 
tion  . 

Rogers  also  investigated  an  enhanced  correlation  tracker 
which  had  "pseudo-measurements"  as  an  output.  The  enhance¬ 
ment  was  due  to  the  current  target  FLIR  image  being  corre¬ 
lated  with  the  intensity  template  described  in  the  previous 
paragraph,  rather  than  with  the  previous  image  data.  The 
"pseudo-measurements"  produced  were  position  offsets  between 
the  target  FLIR  image  and  the  center  of  the  FLIR  FOV  in  each 
of  the  two  FLIR  plane  coordinate  directions.  These  offsets 
were  then  used  as  measurements  provided  to  a  Kalman  filter. 
The  preprocessing  of  the  FLIR  data  by  the  enhanced  corre¬ 
lator  algorithm  allowed  the  use  of  a  linear  Kalman  filter 
instead  of  the  previously  used  EKF .  This  reduced  the  com¬ 
putational  loading  when  compared  to  the  previously  used  EKF. 
Performance  was  generally  comparable  in  terms  of  rms  track¬ 
ing  errors,  with  the  correlator-linear  Kalman  filter  having 
lower  mean  tracking  errors  but  higher  standard  deviations 
than  the  EKF. 

Kozemchak  [3]  and  Millner  [ill  continued  the  algorithm 
development  with  investigations  of  both  the  EKF  and  Roger's 
enhanced  correlator/linear  Kalman  filter,  using  more 
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realistic  target  trajectories.  Both  Gauss-Markov  and 
constant  turn-rate  models  were  considered  for  acceleration 
representation  within  the  filter.  Because  of  the  goal  of 
attempting  to  maintain  tracking  lock  on  harshly  maneuvering 
targets,  the  scheme  used  by  Harnly  and  Jensen  [2]  for 
adaptively  estimating  filter  driving  noise  was  also  imple¬ 
mented.  Both  filtering  schemes  exhibited  good  tracking 
performance  for  target  maneuvers  up  to  five  g's.  However, 
the  abrupt  onset  of  harsh  maneuvers  led  to  serious  perfor¬ 
mance  degradation. 

In  attempting  to  solve  the  problem  of  tracking  harshly 
maneuvering  targets,  Flynn  [1]  proposed  the  use  of  a  multi¬ 
ple  model  adaptive  filter  (MMAF).  Suizu  [17],  following  up 
on  Flynn's  research,  implemented  a  bank  of  two  elemental 
Kalman  filters.  One  filter  was  tuned  for  benign  target 
maneuvers  and  accepted  measurements  from  a  narrow  (8x8 
pixel)  FLIR  FOV.  The  second  filter  was  tuned  for  dynamic 
target  maneuvers  and  used  measurements  from  a  wide  (24  x  24 
pixel)  FLIR  FOV.  Using  a  Bayesian  probabilistic  weighted 
average  of  the  elemental  filter  outputs  [7:129-136],  the 
MMAF  tracker  was  able  to  maintain  lock  on  targets  whose 
dynamics  ranged  from  benign  to  20-g  pull-up  maneuvers  at  20 
kilometers.  The  elemental  filters  were  implemented  using 
both  the  EKF  and  the  enhanced  correlator/linear  Kalman 
filter.  Performance  results  were  comparable,  with  the 
correlator/linear  Kalman  filter  having  smaller  mean  errors 
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and  larger  standard  deviations  than  the  EKF ,  as  seen  in 
earlier  work  of  Rogers. 

Loving  [5]  expanded  Suizu's  HMAF  implementation  by 
adding  a  third  elemental  filter  to  the  MMAF  structure.  This 
filter  was  tuned  for  intermediate  target  dynamics  and  used 
measurements  from  a  narrow  (P  x  8  pixel)  FLIR  FOV.  Loving 
also  compared  a  Maximum  A  Posteriori  (MAP)  MMAF  scheme  to 
the  Bayesian  MMAF  structure  previously  used.  The  MAP  struc¬ 
ture  differs  from  the  Bayesian  structure  in  using  only  the 
estimates  from  the  elemental  filter  with  the  highest  proba¬ 
bilistic  weight,  rather  than  a  weighted  sum  of  elemental 
filter  outputs  as  in  the  Bayesian  MMAF.  The  addition  of  the 
third  filter  enhanced  the  MMAF  tracking  performance.  No 
significant  difference  in  performance  was  found  between  the 
Bayesian  and  MAP  MMAF  algorithms. 

Netzer  [12]  continued  further  research  into  the  three- 
filter  Bayesian  MMAF  algorithm.  After  examining  the  steady 
state  bias  errors  produced  in  the  more  benign  of  the  two 
FLIR  plane  directions  when  the  target  made  a  2U-g  turn,  he 
determined  the  need  to  investigate  an  MMAF  structure  with 
elemental  filters  tuned  for  target  maneuvers  predominantly 
in  either  azimuth  or  elevation.  This  segregation  of  maneu¬ 
ver  direction  into  different  elemental  filters  allows  the 
MMAF  algorithm  to  distinguish  maneuvers  in  these  two  direc¬ 
tions.  This  enables  the  tracker  to  expand  the  FOV  in  the 
critical  direction  on  a  harshly  maneuvering  target.  The 
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tracker  can  then  maintain  lock  on  a  maneuvering  target  in 
the  critical  direction  while  still  producing  accurate  esti¬ 
mates  in  the  direction  of  benign  maneuvers.  Net2er  also 
recommended  the  use  of  a  constant  turn  rate  process  for 
modeling  acceleration  at  close  ranges. 

Tobin  [18]  further  expanded  the  research  into  MMAFs 
using  Netzer's  recommendations .  He  added  two  elemental 
filters,  tuned  for  maneuvers  in  azimuth  and  elevation,  to 
the  three-filter  bank.  A  comparison  between  constant  turn 
rate  and  Gauss-Markov  processes  as  acceleration  models  in 
the  elemental  filters  was  also  made.  Although  the  results 
indicated  smaller  steady  state  errors  for  the  constant  turn- 
rate  model,  the  Gauss-Markov  model  filters  had  consistently 
better  transient  error  performance.  The  performance  analy¬ 
sis  of  the  five-filter-bank  MMAF  indicated  the  tracker's 
ability  to  maintain  target  lock  during  a  jink  in  elevation 
while  still  producing  superior  tracking  performance  in  the 
azimuth  direction  over  an  MMAF  without  any  directionally 
tuned  elemental  filters.  This  was  a  preliminary  feasibility 
study,  with  an  eventual  goal  of  using  an  arbitrary  direction 
for  maneuver  acceleration,  as  opposed  to  strictly  azimuth  or 
elevation . 

Leeney  [4]  continued  the  MMAF  investigation  by  adding 
states  to  the  truth  model  to  account  for  bending/vibration 
effects  in  a  large  space  structure.  Although  bending/vibra¬ 
tion  states  were  not  added  to  the  filter  model,  the 
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elemental  filters  were  retuned  using  the  updated  truth 
model.  The  retuned  MMAF  tracker  exhibited  satisfactory 
tracking  performance  against  targets  undergoing  dynamic  ( 10- 
g)  maneuvers,  provided  that  the  values  given  to  the  para¬ 
meters  associated  with  the  bending/vibration  effects  are  on 
the  order  to  be  expected.  Leeney  also  examined  the  use  of  a 
higher  sampling  rate  (50  Hertz  versus  30  Hertz),  but  found 
the  additional  computational  burden  outweighed  the 
performance  increase.  A  preliminary  investigation  into  the 
use  of  a  rotating  rectangular  (24  x  8  pixels)  field  of  view 
( RRFOV ) ,  in  which  the  elongated  side  of  the  FOV  would  be 
aligned  with  the  estimated  acceleration  direction,  was  done. 
The  rationale  was  to  replace  the  two  rectangular  FOV 
elemental  filters  of  Tobin's  research  with  a  single  filter. 
The  initial  investigation  revealed  that  acceleration  state 
estimates  were  too  noisy  for  accurately  estimating  the 
maneuver  direction;  tne  more  precise  velocity  estimates  were 
used  instead,  with  the  long  side  of  the  rectangular  FOV 
being  aligned  perpendicular  to  the  velocity  vector. 

Leeney's  preliminary  work  did  indicate  enough  promise  with 
this  idea  for  further  exploration 

Norton  [13]  implemented  the  RRFOV  in  his  research.  In 
addition,  he  also  considered  the  effect  of  filter  dynamics 
driving  noise  strength  versus  FOV  size  on  rilter  perfor¬ 
mance.  His  results  showed  that  choice  of  a  larger  dynamics 
driving  noise  strength  in  the  direction  or  maneuver  was  more 
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important  than  increased  FOV  si2e  in  improving  filter  per¬ 
formance.  By  using  an  8  x  8  pixel  rotating  FOV  (instead  of 
an  8  x  24  pixel  rectangular  FOV)  in  combination  with  larger 
values  of  dynamics  driving  noise  in  the  direction  of  maneu¬ 
ver,  an  improvement  in  performance  was  achieved.  Norton 
also  investigated  a  means  to  rotate  the  dynamics  driving 
noise  strength  "Q"  matrix  mathematically  so  that  the  larger 
"Q"  values  stayed  aligned  with  the  acceleration  direction. 

A  scheme  for  simulating  the  physical  rotation  of  the  FLIR 
sensor  plane  to  keep  one  axis  aligned  with  the  acceleration 
vector  was  also  investigated.  Preliminary  results  from 
implementation  of  these  various  schemes  into  an  MMAF  algo¬ 
rithm  were  encouraging. 

Most  recently,  Rizzo  [14]  used  Norton's  results  to 
investigate  a  ballistic  missile  tracking  scenario.  Inter¬ 
ested  in  "pogo"  phenomenon,  where  the  missile  exhaust  plume 
oscillates  along  the  missile  hardbody  longitudinal  axis,  he 
modeled  this  as  a  second  order  process.  Two  parameters, 
natural  frequency  and  oscillation  amplitude,  were  studied, 
and  elemental  filters  from  a  four-filter  bank  were  tuned  for 
different  combinations  of  these  two  parameters.  A  compari¬ 
son  between  a  rotating  FOV  and  a  diagonal  rotating  FOV,  in 
which  the  maneuver  direction  was  aligned  with  the  FOV  diago¬ 
nal  rather  than  one  of  the  FLIR  plane  axes,  was  also  made. 
The  results  indicated  superior  tracking  performance  using 
the  diagonal  rotating  FOV. 
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For  his  research,  Rizzo  used  Norton's  eight-state  ele¬ 
mental  filters  as  benchmarks  for  comparison.  These  filters 
were  then  augmented  to  ten  states  by  the  addition  of  two 
states  associated  with  "pogo"  effects.  Both  the  eight  and 
ten-state  filters  were  then  tuned  against  Norton's  12-state 
truth  model,  which  was  increased  to  14  states  by  the  addi¬ 
tion  of  the  same  2  "pogo"  states.  After  tuning  both  sets  of 
filters,  Rizzo  observed  the  performance  of  the  eight-state 
filters  to  be  superior  to  that  of  the  ten-state  filters. 
Since  the  ten-state  filters  had  knowledge  of  the  "pogo" 
effects  which  were  also  included  in  the  truth  model,  and 
should  have  outperformed  the  eight-state  filters  which  had 
no  knowledge  of  the  "pogo"  effects,  a  preliminary  stochastic 
observability  analysis  [6:243]  was  conducted  to  determine 
possible  causes  for  this  anomaly.  The  analysis  revealed 
that  target  velocity  and  acceleration  states  were  almost 
unobservable  for  the  specific  class  of  missile  targets  he 
was  considering.  Recommendations  for  future  work  included 
reducing  the  order  of  the  filter  models,  and  remodeling  the 
velocity  states  as  first  order  Gauss-Markov  processes. 
Because  of  the  performance  anomalies  discovered,  Rizzo  did 
not  implement  an  MMAF  structure. 

1 . 3  Objectives 

Previous  research  at  AFIT  has  been  directed  toward  the 
goal  of  producing  an  autonomous  tracker.  The  progression  of 
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past  research  shows  a  clear  trend  toward  development  of  an 
algorithm  that  tracks  aircraft  and  ballistic  missiles  under 
increasingly  realistic  conditions.  These  efforts  have 
concentrated  on  tracking  targets  using  only  passively  ac¬ 
quired  infrared  measurements.  Unfortunately,  by  doing  so, 
the  filter  has  no  information  with  which  to  estimate  the 
missile  hardbody  location,  as  separated  from  the  exhaust 
plume  high  intensity  peak.  This  thesis  effort  will  consider 
using  actively  acquired  measurements  as  well,  in  order  to 
help  resolve  the  hardbody  itself. 

Specifically,  the  scenario  investigated  will  consist  of 
tracking  a  ballistic  missile  in  boost  phase  through  the 
atmosphere.  The  ground  sensors  will  consist  of  a  FLIR 
sensor  detecting  exhaust  heat  and  an  additional  optical 
sensor  that  detects  the  returns  from  a  low  power  laser 
illuminating  the  missile  hardbody.  By  observing  the  laser 
returns,  some  additional  information  to  help  discern  the 
missile  hardbody  from  the  exhaust  plume  is  provided  to  the 
filter.  The  specific  objectives  of  this  thesis  effort  are 
outlined  below. 

1.3.1  Observability  Analysis.  Based  on  Rizzo's  recom¬ 
mendations  [I4!/  a  stochastic  observability  analysis  [6:243] 
will  be  conducted.  Because  of  the  observability-like  prob¬ 
lem  encountered  by  Rizzo,  Norton's  eight-state  filter  [13] 
will  be  used  for  this  study.  The  analysis  will  indicate 
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which  states  are  least  observable,  and  thus  the  candidates 
for  elimination  in  an  order  reduction  effort  to  address  this 
issue.  Rizzo  [14:7-4]  indicated  target  acceleration  and 
perhaps  velocity  as  the  two  most  likely  candidates.  Once 
the  appropriate  reduced  order  filter  has  been  determined,  it 
will  be  used  as  the  benchmark  filter  for  this  investigation. 

1.3.2  Active  Illumination  Modelling.  In  the  scenario 
considered  for  this  thesis  effort,  illumination  of  the 
missile  hardbody  is  accomplished  by  sweeping  a  low  power 
laser  up  and  down  the  estimated  missile  velocity  vector, 
starting  from  the  estimated  centroid  of  intensity.  The 
dithering  of  the  laser  would  illuminate  all  or  some  part  of 
the  missile  hardbody.  By  viewing  the  speckle  or  other 
properties  of  the  return  signal,  there  should  be  some  indi¬ 
cation  of  the  line  defining  the  boundary  between  the 
hardbody  and  the  spatial  region  outside  the  hardbody.  For 
this  preliminary  work,  a  modelling  of  the  physical  process 
of  target  illumination  by  a  laser  will  not  be  attempted. 

This  effort  will  instead  concentrate  on  specifying  the  form 
of  the  measurement  data  presented  to  the  filter  from  this 
process . 

The  rationale  for  using  laser  illumination  is  to  provide 
information  to  the  filter  about  the  location  of  the  missile 
center  of  mass  (CM)  relative  to  the  FLIR  plane  intensity 
image  centroid.  The  measurement  form  used  in  this  thesis 
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would  present  this  as  a  measurement  of  the  CM  offset  from 
the  image  centroid  along  the  estimated  velocity  vector.  A 
separate  single  state  filter  will  be  used  to  estimate  the 
offset  distance,  based  on  this  measurement. 

1.3.3  FLIR  Filter.  The  reduced  order  filter  from 
Section  1.3.1  will  be  established  as  the  basic  filter  used 
for  all  investigations  in  this  research.  Before  beginning 
any  of  the  parameter  studies  involving  the  single-state 
filter,  the  reduced  order  filter  will  be  tuned  using  Nor¬ 
ton's  twelve-state  truth  model.  The  truth  model  will  be 
changed  to  reflect  a  ballistic  missile  target,  rather  than 
the  highly  dynamic  aircraft  target  used  by  Norton. 

1-3.4  Single-State  Filter  Performance.  Once  the 
filter  from  Section  1.3.3  has  been  tuned  for  best  tracking 
performance,  a  series  of  parametric  studies  involving  the 
single-state  filter  will  be  performed.  A  sensitivity  study 
will  be  conducted  first  to  determine  the  effect  of  changing 
the  laser  reflection  measurement  noise  level  on  the  single¬ 
state  filter's  performance.  For  this  sensitivity  study,  the 
filter  will  be  correctly  informed  of  the  changes  made  in  the 
"real  world"  tracking  scenario.  Once  the  sensitivity  study 
has  been  completed,  a  robustness  study  involving  reflection 
measurement  noise  level  will  be  performed,  in  which  the 
filter  will  not  be  informed  of  changes  in  the  "real  world". 

A  series  of  parameter  studies,  involving  variation  of  FLIR 
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sensor  sensitivity,  and  separation  distance  of  image  cen¬ 
troid  from  missile  CM,  will  complete  the  analysis. 

1 . 4  Thesis  Overview 

This  chapter  has  reviewed  past  research  and  enumerated 
the  objectives  for  this  thesis  effort.  Chapter  2  considers 
the  truth  model  and  develops  the  simulation  model  used  in 
this  and  past  investigations.  Chapter  3  presents  the  de¬ 
tailed  filter  development.  Chapter  4  develops  the  complete 
algorithm  used  in  this  research  and  discusses  the  manner  in 
which  the  results  of  the  analyses  will  be  presented.  The 
results  of  the  performance  analyses  follow  in  Chapter  5. 

The  final  chapter  summarizes  the  conclusions  and  makes 
recommendations  for  future  study. 


2 .  Truth  Model 


2 . 1  Introduction 

A  truth  model  is  the  result  of  an  attempt  to  describe 
the  behavior  of  real  world  phenomena  of  interest  in  an 
accurate  mathematical  form.  In  many  cases,  the  complete 
description  of  true  system  behavior  may  require  an  infinite 
dimensional  state  space  model.  The  number  of  truth  model 
states  must,  of  necessity,  be  finite  dimensional,  capturing 
enough  of  the  dominant  characteristics  of  system  behavior 
for  an  accurate  representation,  while  at  the  same  time  not 
becoming  too  large  to  be  computationally  unmanageable.  In  a 
simulation,  such  as  the  one  used  in  this  thesis  effort,  the 
truth  model  provides  the  basis  for  depicting  the  true  "real 
world"  dynamic  behavior  of  a  system.  The  truth  model 
becomes  the  standard  against  which  the  performance  of  the 
Kalman  filter  is  evaluated. 

For  this  thesis,  the  true  dynamics  of  the  apparent 
target  image  in  the  FLIR  plane  are  due  to  the  summed  effects 
of  true  target  motion,  atmospheric  jitter,  and  hardware 
vibration/bending.  The  dynamics  are  represented  as  changes 
in  the  FLIR  plane  x  and  y  coordinates  (i.a.,  azimuth  and 
elevation)  of  the  target  intensity  image.  The  position  of 
the  target  image  at  any  one  time  is  given  by: 


=  *t  +  +  xb 


>V  =  Yf  +  +  Yjy 


(  2-1  ) 
l  2-2  ) 
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where 


xc,  y c  =  image  centroid  coordinates 

xt ,  yf  -  component  of  xc ,  yc  due  to  actual  target  dynamics 
xa,  yn  -  component  of  xc,  y due  to  atmospheric  jitter 
x^,  y^  =  component  of  xc,  j /  due  to  mechanical  bending/ 
vibration  of  hardware 

This  chapter  will  describe  the  state  space  structure  of 
this  truth  model.  The  twelve-state  dynamics  model  necessary 
to  account  for  the  terms  in  Equations  (2-1)  and  (2-2)  will 
be  developed  first,  followed  by  the  appropriate  measurement 
model.  The  chapter  will  then  discuss  the  simulation  space 
used  in  this  thesis  effort.  This  discussion  will  include 
descriptions  of  the  coordinate  frames,  target  models,  and 
derivations  of  all  the  measurement  simulations. 

2 . 2  Dynamics  Model 

The  twelve-state  dynamics  model  used  for  the  truth  model 
consists  of  a  two-state  target  dynamics  model,  a  six-state 
atmospheric  jitter  model,  and  a  four-state  mechanical  vibra¬ 
tion/bending  model.  This  model  is  represented  by  the  fol¬ 
lowing  first  order  linear,  stochastic  differential  equation: 

x(t)  =  Fj-Xj- ( l )  +  BjU T{t)  +  w T(t)  (2-3) 

where : 

F-  =  12  x  12  time-invariant  truth  modei  system 

plant  matrix 

B j  =  12  x  2  time  invariant  truth  model 
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input  distribution  matrix 
ytp(t)  =  12-dimensional  truth  model  state  vector 
Ujr( T)  =  2-dimensional  truth  model  input  state  vector 
w j{t)  =  12-dimensional  zero  mean,  white,  Gaussian 

noise  vector  witn  autocorrelation  function: 

E[wT(t)wUt  +  r)  ]=Qt6(t)  .  (2-4) 

The  equivalent  discrete-time  model  [6]  corresponding  to 
Equation  (2-3)  is  of  the  form: 

=  +  B7'du7'rf(M  +  (2-5) 

where  the  state  transition  matrix  $j{t  ,t  t)  is  determined 
by  solving  the  differential  equation  [ 6: 40-41 J: 

=  F T*T{t,ti)  (2-6) 

using  the  initial  condition:  ^(^-,^)  =  I 
The  discrete-time  input  distribution  matrix  in  Equation 
(2-5)  is: 

Brrf  =  t1;1  ^■(tt.rr)&rdT  (2-7) 

i 

Finally,  the  vector  quantities  in  Equation  (2-5)  are: 

x^((,)  =  12-dimensional  discrete-time  truth  model 
state  vector 

u j-tfitj)  =  2-dimensional  discrete-time  input  vector 
Wj-(f(/,)  =  12-dimensional  discrete  time,  zero  mean 

white  Gaussian  noise  vector  with  covariance: 
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lTd 


=  *r(ti-J'T)Q1-tf-{tl.1,T)dT  (2-8) 


where  Qj  is  defined  in  Equation  (2-4 


As  previously  mentioned,  the  truth  model  consists  of  a 
target  dynamics  position  state,  three  atmospheric  jitter 
shaping  filter  [10]  states,  and  two  mechanical  bending 
states,  in  both  x  and  y  FLIR  plane  directions.  In  augmented 
form,  the  truth  model  state  vector  becomes: 


*T 


_  *6  _ 


(  2-9) 


The  discrete  time  truth  model  state  transition  matrix  is: 


*T  - 


*d. 


o 


2  >:  2 \ 


0  2X. 


®(5x2> 


.:-:2  :■ 


(2-10) 


and  the  discrete-time  truth  model  distribution  matrix  is: 


*Td 


B, 


2x2 

0  6x2 

0 


(2-11) 


and  the  discrete-time  truth  model  white  Gaussian  noise 
process  is  given  by: 


w 


Td 


w 


da 


w 


db 


(  2-12  ) 
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where : 


=  2-dimensional  target  dynamics  state  vector 
x^  =  6-dimensional  atmospheric  jitter  state  vector 
x^  =  4-dimensional  bending/vibration  state  vector 


w c/a(^!)  =  6-dimensional  discrete  time,  white  Gaussian 
noise  related  to  atmospheric  states 
w ^(tj)  =  4-dimensional  discrete-time,  white  Gaussian 
noise  related  to  bending  states 

After  examining  Eguations  (2-5)  and  (2-9)  to  (2-12),  the 
truth  model  state  dynamics  representation  is  observed  to  be 
composed  of  three  totally  independent  partitions.  This 
allows  the  dynamics,  jitter,  and  bending  models  to  be  devel¬ 
oped  separately.  The  next  subsections  details  the  develop¬ 
ment  of  these  three  models. 

2.2.1  Target  Dynamic  States.  For  this  and  previous 
thesis  efforts  [1-5,10-18],  target  dynamics  are  modelled  as 
they  occur  in  the  FUR  plane.  The  FLIR  image  plane  is 
modelled  as  being  coincident  with  the  array  of  FLIR  sensor 
elements,  and  perpendicular  to  the  sensor-target  line-of- 
sight  (LOS)  vector.  Since  the  sensor-to-target  range  is 
large,  the  azimuth  and  elevation  displacements  from  the  FLIR 
image  plane  origin  are  directly  proportional  to  "pseudo" 
azimuth  and  elevation  angles.  These  "pseudo"  angles,  alpha' 
and  beta',  as  well  as  the  FLIR  plane,  are  developed  in 
Section  2.3.4. 
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To  convert  the  "pseudo"  angle  measurement  units  of 
microradians  to  the  displacement  distances  in  pixels,  a 
pixel  proportionality  constant  is  used.  This  constant 
represents  the  angular  FOV  of  a  single  picture  element,  or 
pixel.  Continuing  with  Rizzo's  work  [14],  a  value  on  the 
order  of  15  microradians/pixel  was  initially  used  in  this 
thesis  work. 

Using  the  assumptions  that  the  azimuth  and  elevation 
rates  remain  essentially  constant  over  any  sample  period, 
the  discrete  target  dynamics  model  is: 


where : 


(  a  '  )  (  At) 

-  *,«,)  +  1 


vw  -  vt  - 


(  2-13  ) 

(  2-14  ) 


a'(^)  =  c^-- ,  measured  in  micro-radians/second  and 
constant  over  the  time  interval 
dR' 

B'  (t.{)  =  measured  m  micro-radians/second  and 

constant  over  the  time  interval 
=  sample  time  interval  -  t; 
kD  =  pixel  proportionality  constant, 

15  micro-radians/pixel 


Note  the  minus  sign  in  Equation  (2-14).  This  is  due  to  the 
difference  in  the  y  axis  orientations  between  the  FLIR  plane 
coordinate  frame  and  the  inertial  coordinate  frame.  These 
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coordinate  frames  are  shown  in  Figure  2.1  and  developed  in 
Section  2.3.1. 

Arranging  these  equations  in  state  space  form  yields: 


l_r 1  °  “]  r  W  " 

J  L  0  i  J  L  _ 


a'  (l;) 


(2-15)  | 

| 

i 


Using  the  matrices  of  Equation  (2-15)  and  the  block  form 
of  the  overall  truth  model,  the  upper  left  hand  block  of 
Equation  (2-10)  is: 


and  the  upper  block  of  Equation  (2-11)  is: 


*I’d  = 


At 


0 

At 


(2-17) 


and  the  deterministic  input  vector  in  Equation  (2-5)  is: 


uTdl< 1 1 1 


<*'  (f() 
(3  '(t,) 


(  2-18  ) 


The  truth  model  missile  trajectory  simulation  is  a 
continuation  of  that  used  by  Rizzo  [14].  The  underlying 
dynamics  model  [14:3-10]  is  that  of  a  point  mass  influenced 
by  a  thrust  force  and  a  gravitational  force,  with  all  other 
external  forces  assumed  negligible. 
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The  formulation  of  the  truth  model  target  dynamics 
states  in  deterministic  state  space  form  such  as  Equation 
(2-15)  has  some  advantages.  First,  although  the  truth  model 
trajectory  could  have  been  stored  as  FLIR  plane  position 
coordinates  for  each  sample  time  in  the  simulation,  rather 
than  inputting  the  angular  rate  inputs  to  produce  such  a 
position  time  history,  the  present  form  is  more  flexible  in 
supporting  a  wide  variety  of  missile  trajectories.  For 
instance,  if  a  stochastic,  rather  than  deterministic,  dynam¬ 
ics  model  is  desired,  this  can  be  implemented  easily  by  the 
addition  of  a  noise  term  to  Equation  (2-15). 

2.2.2  The  Atmospheric  States.  Using  power  spectral 
density  characteristics,  it  can  be  shown  that  the  atmos¬ 
pheric  jitter  phenomena  can  be  approximated  by  the  output  of 
a  third-order  shaping  filter  driven  by  white  Gaussian  noise 
[10].  The  Laplace  domain  representation  of  the  shaping 
filter  transfer  function  is  given  by  [10:12]: 

*  (s)  KaAB~ 

(2-19) 

(  s+A  )  (  s+B  )  - 

where : 

xa  =  atmospheric  jitter  position  in  one  direction  of  the 
FLIR  image  plane,  the  output  of  the  shaping  filter 
defined  in  Equation  (2-19) 
wa  =  zero  mean,  unit  strength,  white  Gaussian  noise 
a ,7  =  gain,  adjusted  for  desired  jitter  RMS  value 

A  =  break  frequency,  14.14  rad/sec  ’10; 

B  =  break  frequency,  659.5  rad/sec  i 10 < 
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The  inverse  Laplace  transform  of  Equation  (2-19)  is  a 
third-order,  linear  differential  equation.  This  can  be 
expressed  as  three  coupled,  first  order,  linear  differential 
equations  in  state  space  form.  The  atmospheric  jitter  ef¬ 
fects  can  be  modeled  identically  in  both  the  x  and  y  direc¬ 
tions.  The  augmented  six-state  truth  model  for  atmospheric 
jitter  expressed  in  Jordan  canonical  form  is  [10]: 

Xa(0  =  FaX„(n  +  GaWa(t)  (2-20) 

where : 

x^(0  =  6-dimensional  state  vector 

Fa  =  6  x  6  time-invariant  plant  matrix 
wa(/)  =  2-dimensional,  zero  mean,  white  Gaussian 

noise  process  with  unit  strength  and  inde¬ 
pendent  components,  described  as: 

E[wa(n  ]  =  o 

*'[wa(*)wj(*+r)  ]=  Qa6(r)=  l]<5(T) 

The  six  atmospheric  states  in  the  state  vector  correspond  to 
the  low  frequency  pole  and  the  higher  frequency  double  pole 
in  the  x  and  y  FLIR  plane  directions.  The  atmospheric  plant 
matrix  is  defined  in  Jordan  Canonical  form  as  [10:13]: 


~  -A  0  0  0  0  0  " 

0  -B  1  0  0  0 

0  0-5  0  0  U 

0  0  0  -A  U  0 

0  0  0  0  -/?  1 

_0  0  0  U  0  —B  _ 


(2-21) 
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and  the  noise  distribution  matrix  is  [10:13-14]: 


~  Oj  0  " 

G  A  0 


n  G  - 

6 

0  c3 


where : 


(A-B)z 


KnABz 
°3  =  JA-B) 


(  2-22  ) 


The  equivalent  discrete-time  model  associated  with 
Equation  (2-20)  is  the  atmospheric  jitter  partition  of 
Equation  (2-5)  and  is  qiven  by: 

V*-/)  =  +  (2-23) 

Mercier  [10]  showed  that  the  state  transition  matrix  derived 
from  the  time-invariant  plant  matrix  of  Equation  (2-20)  is: 


(A')  = 


♦a//  °  ° 

U  ^<72(2  ®aAJ 


0 

0 

0 

0 


0 

0 

0 

0 


aJJ 

0  0 

0 
0 


0 
0 
0 

all 
0  0 


0  0 
0  0 
0  0 
0  0 


0 


(20  0 


(  2-24) 
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where : 


*all 

*a44 

=  exp(-AAt) 

*<122 

*a55 

=  exp(-BAt) 

*1 23 

*a56 

-  Atexp(-BAt) 

*a33 

*066 

=  exp(-BAt) 

The  six-dimensional,  zero  mean,  discrete-time,  white, 
Gaussian  noise  w  has  statistics  defined  as: 


=  0 


( 2-25a) 


=  Qda  =  J*J+/  *a(‘M'*)<W£*fa(tM,T)dT  (2-25b) 

Li 


2.2.3  Bendinq/Vibration  States.  Leeney  [4]  recently 
added  mechanical  bending  states  to  the  truth  model.  These 
states  were  added  to  account  for  vibrational  effects  in  the 
FLIR  data  that  occur  when  the  sensor  is  mounted  on  a  moving 
non-rigid  platform.  Based  on  AFWL-conducted  tests,  Leeney 
concluded  that  bending  effects  in  both  the  x  and  y  FLIR 
directions  could  be  represented  by  a  second  order  shaping 
filter,  driven  by  white  Gaussian  noise.  The  Laplace  domain 
transfer  function  of  this  shaping  filter  is  represented  by 
[5:33]: 

*b[s_\  =  Kh(i)nb  _  to-of 

Wb(s)  si  +  2hU}nbs+*nb 

where : 

=  shaping  filter  output,  the  FLIR  plane  positional 
offset  due  to  bending/vibration 
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wb  =  zero  mean,  unit  strength,  white  Gaussian  noise 
with  an  autocorrelation  of: 

£['"b(t)wi(t  +  T)  }  =  Qb6(t-T);  Qb  =  1 
Kb  =  gain  adjustment  to  achieve  desired  root  mean 

sguare  (RMS)  bending  output;  Kb  =  5x10  ' 

(Note:  Kb  is  given  here  because  the  strength  of  the 
bending  white  noise  is  expressed  in  terms  of  this 
parameter,  rather  than  Kb) 

=  damping  coefficient,  =  0.15 
unb  =  undamped  natural  bending  freguency,  n  rad/sec 

Leeney  [4:35]  determined  that  the  x  and  y  directions  could 
be  treated  independently  and  augmented  them  together  to  form 
a  four-state  model.  The  linear  stochastic  differential 
eguation  describing  the  bending/vibration  effects  is: 

x&(0  =  F6V')  +  (2-27) 

where : 

x^(n  =  4-dimensional  state  vector 

=  4  x  4  time-invariant  plant  matrix 
w^(0  =  2-dimensional,  zero  mean,  white  Gaussian 

noise  process  with  independent  components  of 

r 1  0 

unit  strength,  Qh  = 

L  0  1 

=  4  x  2  noise  distribution  matrix 
The  bending  plant  matrix  is  defined  as  [4:142]: 
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0  1 

~wnb  ~2hwnb 
0  0 

0  0 


0  0 
0  0 
0  1 


-u) 


nb 


-2Cfcu> 


b  nb 


(  2-28  ) 


and  the  noise  distribution  matrix  is  [4:142]: 


05  = 


0 

*nb2kp 


0 

0 


0 

0 


*nb  kp 


(  2-29  ) 


(Note  that  kp  is  the  pixel  proportionality  constant) 

The  equivalent  discrete-time  model  for  Equation  (2-27)  is 
of  the  form: 

=  +  Wdft(<l.)  (2-30) 

where : 


%(&) 


and 


"  *&7 

*b2 

0 

0  " 

*5J 

hi 

0 

0 

n 

0 

hi 

<t>b2 

0 

0 

h3 

hi  . 

<t>bJ  =  exp ( -Opht )  [  cos  ( up&t )  +  ^-sin(w6W)] 
<t>b2  =  exp(-aftAO  [^sin(wftW  )  ] 

V;  =  exp(-a6AH  [-l-(^)^  ]sin(w6A? ) 

a6  ■ 

<t>b4  =  exp (  -o^t )  [  cos  ( ijpM  )  -  y-  sin(w44n] 
A7  =  sample  time  interval  [  f. ,  -  /,■] 


(2-31) 
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a *  =  real  part  of  the  root  of  the  characteristic 
equation  in  Equation  (2-26) 

=  imaginary  part  of  the  root  of  the 

characteristic  equation  in  Equation  (2-26) 

The  four-dimensional,  zero  mean,  white  Gaussian  noise  pro¬ 
cess  of  Equation  (2-30)  has  an  equivalent,  discrete-time 
representation  with  statistics  given  by: 


f[wd5(^)]  =  0 


(  2-32 


)»'„((,)!  -  Qrfi,  -  fj'1  VWT>WW<(1-/'r>‘ir  12-33) 


TJ, 


2 . 3  Simulation  Space 

In  order  to  simulate  the  operation  of  a  FLIR  sensor 
accurately,  a  "simulation  space"  has  been  developed  for  use 
on  a  digital  computer.  This  simulation  space  has  two 
purposes.  First,  the  representation  of  a  realistic  target 
trajectory  through  three-dimensional  space  can  be  achieved. 
Second,  the  simulation  space  provides  the  means  of 
mathematically  translating  the  target  infrared  image  and 
velocity  vector  in  three-dimensional  space  onto  the  two- 
dimensional  FLIR  image  plane.  These  translations  will  be 
discussed  in  this  section.  However,  the  coordinate  frames 
which  provide  the  basis  for  these  transformations  will  be 
presented  first. 


2.3.1  Coordinate  Frames.  The  following  coordinate 
frames  are  used  in  the  simulation  of  the  operation  of  the 
FLIR  sensor  on  a  digital  computer  [ 14  :  3  :  23-24  ]  : 
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Inertial  Frame: 


Origin:  location  of  the  FLIR  sensor 

Axes:  e^.  -  due  north,  tangent  to  the  earth's  surface, 

defines  zero  azimuth 

&y  -  inertial  "up"  with  respect  to  flat  earth 
approximation 

-  vector  completing  right-hand  coordinate  set, 
defines  90  azimuth 

Note:  The  azimuth  angle  (a)  is  measured  eastward  from  ex. 

The  elevation  angle  (6)  is  measured  "up"  from  the 
horizontal  plane  defined  by  ex  and  . 

Target  Frame: 

Origin:  center  of  mass  of  the  target 

Axes:  ey  -  along  the  true  velocity  vector 

epu  -  out  the  right  side  of  the  target, 
perpendicular  to  ey 

eppu  ”  v^ctor  completing  the  right-hand 
coordinate  set 

Note:  ' u'  -  along  the  velocity  vector 

' p u'  -  perpendicular  to  the  velocity  vector 

'  ppu'  -  perpendicular  to  both  of  the  above 

a  -  6  -  r  Frame: 

Origin:  center  of  mass  of  the  target 

Axes:  ey  -  coincident  with  the  true  sensor-to-target 

LOS  vector. 
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ea,  -  define  a  plane  perpendicular  to  er, 

rotated  from  inertial  and  by  the 
azimuth  angle  (a)  and  elevation  angle  ( B)  . 

FLIR  Image  Plane: 

This  is  the  FLIR  image  plane  defined  by  the  eQ  and 
unit  vectors.  Because  of  small  angle  approximations, 

the  "pseudo"  azimuth  and  elevation  angles,  a' and  3'  as 
shown  in  Figure  2.1,  measured  with  respect  to  the  FLIR  LOS 

vector,  are  linearly  proportional  to  the  x  and  y  cartesian 
coordinates  on  the  FLIR  plane.  The  x  and  y  coordinates  are 
distances,  in  pixels,  from  the  center  of  the  FLIR  FOV. 
Observing  the  FLIR  plane  from  the  inertial  origin,  x  is 
positive  to  the  right  and  y  is  positive  down.  This  conven¬ 
tion  is  used  to  maintain  a  right-handed  coordinate  system. 

The  inertial  and  target  frames,  as  well  as  the  FLIR  image 
plane,  are  illustrated  in  Figures  2.1.  and  2.2.  Note  that, 
for  interpreting  the  projections  in  Figure  2.2,  the  e, 
direction  is  not  necessarily  in  the  FLIR  image  plane;  this 
will  be  developed  further  in  Section  2.3.3. 

2.3.2  Target  Model  The  basic  target  model  used  for 
this  and  previous  theses  [14]  is  a  planform  with  two  inten¬ 
sity  functions.  The  spatial  relationship  between  the  two 
intensity  functions  is  shown  in  Figure  2.3.  The  displace¬ 
ment  of  the  two  Gaussian  intensity  function  centroids  along 
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Figure  2  1  The  FLIP  Plane 
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the  e,  axis  was  based  on  rhe  assumption  that  the  dispersion 
of  the  exhaust  plume  in  the  e?v  direction  is  approximately  20 
times  the  radius  of  the  missile  [14].  The  centroid  of  the 
first  intensity  function  is  located  b5  meters  behind  the 
center  of  mass  of  the  missile.  This  distance  was  chosen  tc 
simulate  the  composite  centroid  of  the  exhaust  plume  being 
close  to  the  missile  exhaust  nozzle,  using  the  assumption 
that  the  distance  from  the  missile  center  of  mass  to  the  end 
of  the  missile  is  20  meters.  The  second  intensity  function 
is  located  110  meters  from  the  center  of  mass.  The  rela¬ 
tionship  between  the  missile  center  of  mass  and  the  centers 
of  these  intensity  functions  remain  fixed  in  the  target 
frame  during  a  simulation,  and  is  indicated  in  Figure  2.3. 

As  was  previously  indicated  in  Section  2.2.1,  any  external 
forces  acting  on  the  missile  other  than  thrust  and  gravita¬ 
tional  forces  are  assumed  negligible.  Therefore,  sideslip 
angle  and  angle  of  attack  are  considered  to  be  zero.  These 
simplifications  allow  the  semi-major  axes  of  the  elliptical 
constant-intensity  contours  of  each  of  the  infrared  inten¬ 
sity  functions  to  be  aligned  with  the  target's  velocity 
vector,  as  illustrated  in  Figure  2.4.  As  noted  in  previous 
thesis  efforts  [12],  this  yields  a  simplification  of  the 
simulation  space  geometry  while  retaining  the  essential 
features  of  the  trajectory  simulation  necessary  for  the 
performance  analysis  of  the  tracker. 
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2.3.3. 


During  the  simulation,  as  the  target  propagates  through 
three-dimensional  inertial  space,  the  output  of  the  FLIR 
detector  elements  in  the  array  is  simulated  by  projecting 
the  target's  two  intensity  functions  onto  the  FLIR  plane. 

In  this  and  previous  research  [13,14],  each  of  the  intensity 
functions  remains  fixed  with  respect  to  the  target  frame, 
while  the  locations  of  the  intensity  functions  in  the  FLIR 
frame  change  as  the  target's  orientation  relative  to  the 
sensor  also  changes.  For  simplicity,  the  location  of  each 
of  the  intensity  functions  (hotspots)  is  initialized  in  the 
target  frame  as  a  displacement  from  the  missile  center  of 
mass  (Figure  2.3).  To  orient  the  intensity  functions  in  the 
FLIR  coordinate  frame,  they  are  rotated  by  the  target  orien¬ 
tation  angle  &j-  (Figure  2.4). 

Consider  the  geometry  of  Figure  2.5.  This  figure  shows 
the  geometrical  relationship  between  the  current  target 
image  and  a  "reference  target"  image  in  the  FLIR  plane 
(Figure  2.3  or  2.4).  The  reference  image  is  oriented  to 
correspond  to  the  largest  apparent  planform  at  a  given 
range.  Using  Figures  2.4  and  2.5,  the  current  single  hot¬ 
spot  image  is  defined  as  [5:37]  : 

apu  =  apuolV  )  (2-34) 
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where : 

au0,  apu0  =  the  initial  dispersion  of  the  target 

intensity  function  along  ^  and  epu  in  the 

target  frame  of  the  reference  image 
°u 1  °pv  ~  the  current  dispersions  of  the  target  image 

r-  -  initial  sensor-to-target  range  of  the 
reference  image 

r  =  current  sensor-to-target  range 
v  =  target  inertial  velocity  vector 
u  =  magnitude  of  v 

v  =  projection  of  v  onto  the  FLIR  plane; 

the  component  of  v  perpendicular  to  the 
LOS  vector 

v .  10$  =  magnitude  of  v^-;  u  /^s=\a"  +  fi 
y  =  angle  between  v  and  the  FLIR  plane 

AR  =  :  aspect  ratio  of  the  reference  image 

puo 

Eguations  (2-34)  and  (2-35)  define  the  dispersion  along  the 
principle  axes  of  the  intensity  function's  constant-inten¬ 
sity  ellipses  as  seen  by  the  FLIR  sensor  (Figure  2.4).  Such 
a  description  is  accomplished  for  each  of  the  two  intensity 
functions  of  Figure  2.3. 

2.3.4.  Velocity  Projection  onto  the  FLIR  Plane . 

The  deterministic  input  vector,  u =  h'(tl)\t  in 

Equation  (2-5),  is  the  projection  of  the  target's  inertial 
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velocity  vector  onto  the  FLIR  image  plane.  Harnly  and 
Jensen  [2]  demonstrated  that  this  projection  is  based  on  the 
geometry  illustrated  in  Figure  2.6.  From  Figure  2.6,  it  can 
be  seen  that: 


a(t)  =  arctan  £  (2-36) 

Taking  the  time  derivative  of  Equation  (2-36)  and  using 
the  fact  that  the  sensor-to-target  range  is  large  so  that 

a(t)  =  a'  (t)  yields: 


a  '  ( t )  =  a ( t) 

where : 


x(t)vz(t)- z(t)vx ( t ) 
xz(t)+zz(t) 


(  2-37  ) 


uA,,uz  =  components  of  the  target's  inertial  velocity 
in  the  ex  and  ez  directions 


In  a  similar  development: 


6(0  =  arctan 


B'(()  =  6(0  = 


rnni 
L  vo  J 


VO 

rh{t)\>y{t)-y{t)'rh(t) 

0(0 


(  2-38 


(  2-39 


where : 


0(0  = 


^(0u,(0^(0u7(0 

VO 


u y  =  component  of  the  target's  inertial  velocity  in 
the  ey  direction 


Equations  (2-37)  and  (2-39)  define  the  deterministic  input 
vector  u r-yO,)  in  the  truth  model  dynamics  difference 
equation,  Equation  (2-6). 
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2 . 4  Measurement  Models 


The  realistic  modelling  of  measurements  of  real-world 
phenomena  is  an  important  feature  of  any  simulation.  This 
section  of  the  chapter  will  describe  the  method  for  simulat¬ 
ing  the  measurement  output  from  each  infrared  detector 
element  in  the  FLIR  array  which  has  been  used  in  previous 
efforts,  as  well  as  in  this  effort.  The  method  for  simulat¬ 
ing  the  output  from  an  optical  sensor  which  receives  low 
power  laser  light  reflected  from  the  missile,  as  developed 
for  this  thesis  effort,  will  also  be  discussed. 

2.4.1  Infrared  Measurement  Model  Information  on  the 
target  is  obtained  by  measuring  the  average  intensity  of  the 
infrared  energy  received  by  each  detector  element  in  the 
FLIR  array.  The  target's  apparent  infrared  image  or  "inten¬ 
sity  function"  on  the  pixel  array  of  detectors  is  the  col¬ 
lective  sum  of  effects  due  to  target  exhaust  plume  infrared 
radiation,  background  noise,  and  sensor  noise. 

Consider  the  energy  radiated  from  a  target  with  a  single 
intensity  function.  The  infrared  intensity  function  on  the 
FLIR  image  plane  can  be  modelled  as  a  bivariate  Gaussian 
distribution  with  elliptical  constant  intensity  contours 
[2].  This  bivariate  Gaussian  intensity  function  is  given  by 
the  following  eguation  [13]: 

7f  K’V ’Kpuak{t) ’v  peak{t)]  =  7«a.*exp  <  -0 . 5  [  AxAk  ]  P  ‘  [  A  xtiv  }  7  !  (2-40) 
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where : 


=  ^“V)cos0^  +  ^“Va*)sin0r 

Ay  =  ^"V^)COS07'  ~  ^'xpeak^sin&r 

Qj ■  =  target  orientation  angle  between  the  projection 
of  the  velocity  vector  onto  the  FLIR  plane 
and  the  FLIR  plane  *  axis;  See  Figures  2.5  &  2.7 
x,y  =  reference  coordinate  axes  on  the  FLIR  plane 

xpeak '  V peak  =  coordinates  of  the  peak  intensity  of  the 

single  Gaussian  intensity  function 
1  max  =  maximum  intensity  of  the  function 

P  =  2  x  2  target  dispersion  matrix  whose 

eigenvalues  ( and  apu)  define  the  dispersion 

of  the  elliptical  constant  intensity  contours 
(along  the  velocity  vector  and  perpendicular 
to  the  velocity  vector)  in  the  FLIR  plane 
(see  Section  2.3.3) 

The  composite  FLIR  plane  image  intensity  function, 
represented  as  the  difference  between  two  individual  inten¬ 
sity  functions  for  a  missile  exhaust  plume,  is  shown  in 
Figure  2.7.  To  form  the  characteristic  crescent  shape  of  a 
missile  plume,  the  rear  individual  intensity  function  is 
subtracted  from  the  forward  intensity  function.  Since  the 
intensity  value  from  an  FLIR  sensor  element  cannot  be  nega- 
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Figure  2  7  FLIR  Image  Plane  Intensity  Function  for  the 
Difference  Between  Two  Gaussian  Intensity  Functions 
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tive,  the  simulation  software  sets  any  calculated  negative 
intensity  values  to  zero. 

The  intensity  measurement  produced  by  each  sensor  pixel 
is  the  average  infrared  energy  intensity  on  that  pixel. 

This  intensity  is  the  sum  of  the  target's  intensity  func¬ 
tion,  spatially  correlated  background  noise,  and  FLIR  sensor 

noise.  The  output  of  the  sensor  pixel  in  the  row 

and  kth  column  of  the  array  at  time  t,  is  [14  1: 


zjk^0  =  {  peak. }U^yp,akM0  1 

xpeakJti^  '  Ypeak^1^  j  >  d*  dy 


njk^0  +  hlk^i) 


(  2-41 


where : 

=  output  of  pixel  jk ,  average  intensity 
on  that  pixel 
Ap  =  area  of  one  pixel 

/y,/?  =  intensity  function  of  the  first  and  second 


x,y 

xpeak  '  y ppqk 
xpeak  '  ^ peak 


Gaussian  intensity  function  (see  Figures  2.3 
and  2.7) 

coordinates  of  any  point  within  pixel  jk 
coordinates  of  maximum  intensity  of  the  first 

Gaussian  intensity  point 

coordinates  of  maximum  intensity  of  the  second 
Gaussian  intensity  point 

effect  of  internal  FLIR  sensor  noise  on  pixel  jk 
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bjklti'  ~  eff  ect  of  spatially  correlated  background  noise 
on  pixel  jk 

The  sensor  error  n is  the  result  of  thermal  noise 
and  dark  current  in  the  infrared  detectors.  This  error  is 
assumed  to  be  both  temporally  and  spatially  uncorrelated 
[14]  . 

The  background  noise  h;j<  (l t  )  is  represented  as  a 
spatially  correlated  noise  with  radial  symmetry,  whose 
correlation  decays  exponentially.  Harnly  and  Jensen  [2] 
used  a  correlation  distance  of  approximately  twro  pixels  in 
the  FLIR  plane,  and  simulated  this  by  maintaining  non-zero 
correlation  coefficients  between  each  pixel  and  its  nearest 
two  neighbors  in  all  directions. 

By  concatenating  all  64  values  of  6 (corresponding  to 
an  8  x  8  pixel  field  of  view)  into  a  64-dimensional  vector 
b(  t- )  ,  the  spatially  correlated  background  noise  is 
modelled  as  [ 14 ] : 

b(f,-)  =  NR  b'  [t{)  (2-42) 

where : 

R  -  64  x  64  correlation  matrix  of  the  discrete, 

zero  mean,  white  Gaussian  vector  noise  process  b( t ) 
b'(*,)  =  64-dimensional,  discrete,  zero-mean,  white 
Gaussian  vector  noise  process  with  the 
correlation  matrix  I 

=  Cholesky  square  root 


A  detailed  development  of  this  spatially  correlated 
noise  process  as  well  as  the  FLIR  sensor  noise  process  can 
be  found  in  the  work  of  Harnly  and  Jensen  [2J.  It  is 
mentioned  here  only  for  completeness  in  describing  the  truth 
model.  The  actual  parameters  used  for  the  sensor  error  and 
the  background  noise  will  be  detailed  in  Chapter  4. 

2.4.2  Laser  Reflection  Measurement  Model  For  this 
initial  effort,  no  attempt  was  made  to  model  the  physical 
phenomena  involved  with  laser  light  reflected  from  the 
missile  hardbody  being  received  by  an  optical  sensor.  Since 
the  purpose  of  reflected  laser  light  would  be  to  provide 
information  on  the  location  of  the  missile  center  of  mass, 
an  attempt  was  made  to  model  the  information  about  the 
center  of  mass  which  would  be  derived  from  reflected  laser 
light . 

The  scenario  for  development  of  the  missile  center  of 
mass  measurement  is  a  continuation  of  that  used  by  Rizzo 
[14].  A  ballistic  missile  in  boost  phase  is  tracked  using 
measurements  from  the  FLIR  sensor.  The  "pseudo"  measure¬ 
ments  derived  from  the  FLIR  sensor  measurements  and  enhanced 
correlator  are  input  to  a  Kalman  filter,  which  provides 
position  and  velocity  estimates  of  the  infrared  image  cent¬ 
roid  of  the  exhaust  plume.  The  enhanced  correlator  will  be 
described  in  Chapter  4.  These  estimates  are  then  used  to 
aim  a  low  power  laser  along  the  filter  estimated  velocity 
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vector  from  the  estimated  center  of  intensity.  Any  reflec¬ 
tions  received  from  the  missile  provide  information  about 
the  missile  center  of  mass.  By  observing  the  characteris¬ 
tics  of  the  return  signal  from  the  laser  being  swept  along 
the  velocity  vector,  one  can  anticipate  getting  an  indica¬ 
tion  of  when  the  laser  passes  the  interface  between  the 
hardbody  and  just  background,  at  both  ends  of  the  missile 
hardbody.  Since  the  image  intensity  centroid  position  and 
velocity  estimates  are  in  terms  of  FLIR  plane  variables,  a 
decision  was  made  to  simulate  the  mass  center  measurement  in 
the  FLIR  plane  also.  With  this  decision,  the  problem  of 
simulation  then  becomes  a  problem  involving  plane  geometry. 

The  first  phase  of  the  simulation  involves  the  projec¬ 
tion  of  the  three-dimensional  missile  hardbody  onto  the  two- 
dimensional  FLIR  plane.  For  simplification,  a  rectangle  was 
chosen  to  represent  the  shape  of  the  missile  hardbody  pro¬ 
jection.  With  the  assumption  that  the  longitudinal  axis  of 
the  missile  hardbody  is  aligned  with  the  velocity  vector, 
the  geometry  for  projection  of  the  missile's  three-dimen¬ 
sional  length  onto  the  FLIR  plane  is  shown  in  Figure  2.5. 

The  geometry  is  described  by: 

Lengthf-^  =  cosy  Length^^  (2-43) 

where : 

Lengthy ^  =  FLIR  plane  projection  of  missile  length 

y  =  angle  between  v  and  the  FLIR  plane  (radians) 
Lengthy  t)]  .  =  true  missile  length  in  pixels 
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Since  the  missile  is  cylindrical,  the  projection  of  the 
missile  diameter  onto  the  FLIR  plane  is  equal  to  the  missile 
diameter.  Once  the  projection  of  the  missile's  dimensions 
onto  the  FLIR  plane  is  accomplished,  the  rectangle  repre¬ 
senting  the  missile  is  located  on  the  FLIR  plane  by  offset¬ 
ting  the  rectangle's  center  from  the  truth  model  image 
centroid  along  the  FLIR  plane  truth  model  velocity  vector. 
The  offset  distance  was  chosen  by  using  the  intensity  func¬ 
tion  displacements  in  Figure  2.3. 

The  next  step  in  the  simulation  involves  the  determina¬ 
tion  of  the  missile  center  of  mass  measurement.  This  is 
accomplished  geometrically.  The  path  along  which  the  laser 
would  be  dithered  is  represented  by  a  rectangle  (a  rectangle 
rather  than  a  line  was  used  to  account  for  the  laser  beam 
width).  If  the  laser  hits  the  missile,  some  form  of  reflec¬ 
tion  is  obtained.  This  reflection  is  geometrically  simula¬ 
ted  by  having  the  two  rectangles  overlap  as  shown  in  Figure 
2.8a.  The  coordinates  of  this  overlap  area  are  determined  by 
finding  the  intersections  of  the  line  segments  representing 
the  sides  of  the  rectangles.  These  intersection  points,  as 
shown  in  Figure  2.8b,  (representing  the  corners  of  some 
enclosed  area)  are  used  to  determine  the  centroid  of  the 
area  of  intersection  (and  presumably  reflection).  Once  the 
center  of  mass  measurements  coordinates  have  been  calcu¬ 
lated,  the  offset  distance  from  the  filter  estimated  image 
centroid  is  determined.  In  the  simulation,  if  the  combina- 
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tion  of  estimated  intensity  centroid  position  and  estimated 
velocity  vector  orientation  is  not  accurate  enough  to  cause 
the  two  rectangles  to  overlap,  then  no  "reflection"  occurs, 
and  no  measurement  update  occurs  in  the  filter.  In  this 
case,  the  filter  simply  propagates  the  estimate  until  a 
reflection  does  occur.  A  complete  development  of  the  geome¬ 
trical  simulation  of  the  mass  center  measurement  is  detailed 
in  Appendix  A. 

The  offset  measurement  which  is  available  for  filter  use 
is  modeled  by  adding  measurement  noise.  The  output  of  the 
sensor  at  time  t-  is: 

=  L  off  set  +  v«i)  (2- 44) 

where : 

z(t-)  =  sensor  output  at  time  (i-) 

L offset^1 0  ~  offset  of  center  of  mass  from  filter  estimated 
FLIR  image  centroid  along  the  filter-estimated 
velocity  direction  in  the  FLIR  plane 
i /((  )  =  discrete,  zero  mean,  white  Gaussian 

measurement  noise  with  statistics  t  (  i/(  t , )  \  =0  , 

F  i  v(  f  )  v{  l  ( )  ;  =R(  tt )  ,  i=i,  o  otherwise 

The  simulated  measurement  developed  in  this  section  will  be 
related  to  the  linear  filter  measurement  model  in  Chapter  3. 
I'he  parameters  used  tor  the  measurement  noise  strength  will 
be  discussed  in  chapter  3. 
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This  chapter  has  developed  the  truth  model  dynamic 
system  as  the  augmentation  of  a  deterministic  target  trajec¬ 
tory  component,  a  stochastic  component  due  to  atmospheric 
jitter,  and  a  stochastic  component  due  to  mechanical  bend¬ 
ing/vibration.  In  order  to  simulate  the  operation  of  a 
tracker  on  a  computer,  a  simulation  space  with  various 
coordinate  frames  and  a  target  model  was  developed. 

Finally,  a  measurement  model  for  the  output  of  the  infrared 
detectors  in  the  FLIR  array  and  a  model  for  the  output  of 
the  optical  sensor  receiving  returns  from  the  reflected 
laser  light  were  developed. 
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3 .  Filter  Models 


3 . 1  Introduction 

This  chapter  presents  the  filter  models  used  in  this 
thesis  effort.  Section  3.2  develops  the  six-state  filter, 
based  on  FLIR  measurements,  which  consists  of  target  dynamic 
states  and  atmospheric  jitter  states.  Section  3.3  develops 
the  single  state  filter,  based  on  the  processing  of  low- 
power  laser  returns,  which  estimates  the  offset  distance 
between  the  image  centroid  and  the  missile  center  of  mass. 
The  measurement  models  used  in  each  of  the  filters  are  also 
developed . 

3.2  FLIR  Filter 

The  filter  initially  used  at  the  beginning  of  this 
thesis  was  based  on  the  eight-state  benchmark  filter  used  by 
Rizzo  [14].  Based  on  the  results  of  the  observability 
analyses  (see  Chapter  5),  the  target  acceleration  states 
were  deleted  from  the  filter.  For  the  remainder  of  this 
thesis  work,  a  six-state  linear  Kalman  filter  was  used  to 
provide  state  estimates  for  the  FLIR  image  position  and 
velocity,  as  well  as  estimates  of  atmospheric  jitter. 

3.2.1  Dynamics  Model.  Previous  AFIT  research  has 
considered  two  different  models  for  representing  target 
dynamics.  The  first  model  represents  target  acceleration  as 
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a  zero  mean,  first-order  Gauss-Markov  process;  the  second 
describes  acceleration  using  a  constant  turn-rate  model. 

With  the  deletion  of  acceleration  states  from  the  previously 
used  eight-state  filter,  a  new  representation  for  target 
velocity  was  needed.  A  decision  was  made  [8,14]  to  repre¬ 
sent  target  velocity  as  a  zero  mean,  first-order  Gauss- 
Markov  process.  This  is  because  time-correlated  physical 
states  such  as  position,  velocity,  and  acceleration  can 
often  be  well  represented  by  an  exponentially  time-corre¬ 
lated  first-order  Gauss-Markov  process.  The  state  vector 
for  the  six-state  Kalman  filter  used  in  this  research  is 
defined  as: 


where  the  states  are : 

=  x  (i.e.,  azimuth)  component  of  FLIR  image  centroid 

position  due  to  target  dynamics 
y !  =  y  (i.e.,  elevation)  component  of  FLIR  image 

centroid  position  due  to  target  dynamics 
vx  =  x  component  of  FLIR  image  velocity 

vy  -  y  component  of  FLIR  image  velocity 

xn  =  x  component  of  atmospheric  jitter 

y  =  v  component  of  atmospheric  litter 
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Each  element  in  Equation  (3-1)  is  coo rdinatized  in  the 
FLIR  plane.  Note  that  the  atmospheric  jitter  model  has  been 
reduced  from  the  six-state  form  in  the  truth  model  (see 
Section  2.2.2)  to  the  two-state  model  seen  here.  The  effect 
of  the  higher  frequency  double  pole  was  negligible  [13]  and 
was  disregarded  to  reduce  the  filter  order.  Also  note  the 
omission  of  any  bending/vibration  states.  Leeney  [4]  tested 
the  omission  of  bending  states  from  a  similar  filter,  and 
found  no  significant  degradation  in  filter  performance.  The 
filter  model  is  described  by  the  following  time-invariant, 
linear  stochastic  differential  equation: 


Xf(  t )  =  FfXf(t)+GfV  f(t)  (3-2) 

where : 

Xj(0  =  six-state  filter  state  vector  of  Equation  (3-1) 

Ff  =  time-invariant  system  plant  matrix 
G j-  -  time-invariant  noise  distribution  matrix 
w f{t)  =  zero  mean,  white  Gaussian  noise  vector  of 
strength  Q ^ 


and 
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where : 

Tx>  ry  ~  correlation  times  for  the  target  x  and  y 
velocities 

ra  =  correlation  time  for  the  atmospheric  jitter 
position  process 

crx,  o~y  =  variance  and  mean-squared  value  for 
the  target  x  and  y  velocities 
cry  -  variance  and  mean-squared  value  for  the 
atmospheric  jitter  position  process 


The  parameters  used  for  the  filter  correlation  times  and 
variances  will  be  described  in  Chapter  5. 

The  filter  state  estimate  and  error  covariance  matrix 
are  propagated  forward  over  a  sample  period  using  the 
following  equations  16:171-172]: 
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x7(f;w)  =  (3-6) 

p/(  #;-./)  =  *rf  A/ )P/(#~-)^(  Af  )  +  Qrf/  (3-7) 

where : 


Xy(^)  =  filter  estimate  of  6-dimensional  state  vector 
P/(t ,)  =  filter  covariance  matrix  (6x6) 

(t,)  =  time  instant  before  FLIR  measurement  is  incorporated 
into  the  estimate  at  time  t- 

{t~A  =  time  instant  after  FLIR  measurement  is  incorporated 
into  the  estimate  at  time  t- 
=  time  invariant  state  transition  matrix  for 
propagation  over  the  sample  period:  -  t- 

and  Qdf  is  defined  as: 


=  j*;*' 

i 

with 


*/(  AO 


where : 
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(  3-9) 


0;/  =  Tx(  l-exp(— T  )  ) 
=  Tyd-exp(-^)  ) 
0,,  =  exp(-f-) 
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and 


*H  =  exp(-") 

*iS  =  *66  =  exP<-f'> 


q  dfli  0  qdfl3  0  0  0 

0  qdfZZ  0  qdfZ3  0  0 

qdf3!  0  qdf33  00  0 

QdJ  0  qdf32  0  qdf-13  0  0 

oooo  ?rf/55  0 

ooooo  q.t,, 

L  dfoo  J 

where : 

qdfn  =  2o^ tx[M-2tx(  1-exp ( )  )+-^(  1-exp (-2Td/  ))  J 

T 

=  2dj,rJ/(«-2r)/(l-exp(  -f'  )  )+  ■£{  l-exp( -2“  )  )  ] 
=  2<4[2T^(l-exp(-f'))  -  T|(l-exp(-^))] 

T 

qdfZ3  =  2op[  2Ty(  l-exp(  ~jy)  )  -  2y(l-exp(-2^  )  )  j 

qdfV  =  qdf?l 
q  df  33  =  ^U-exp(-2^)) 
qdf32 
q  df  33 

a.,,  =  q  ,r.  =  cr, ( 1  -exp (-%“")  ) 

■ uj /n  y  dfoo  a  c  'a 


qdf73 

.  .  ,  2l\t  ,  , 

=  crr,(  l-exp(  -  T  )) 
k  v 


(  3-10) 


The  pointing  controller  used  for  this  and  previous 
studies  is  idealized.  Physical  implementation  considera¬ 
tions,  such  as  servo  lag  and  inertia,  are  neglected.  Netzer 
;  12;  demonstrated  that  any  errors  generated  by  Li.e  use  of  an 
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idealized  controller  are  small  and  are  interpreted  by  the 
Kalman  filter  as  atmospheric  jitter.  Following  the  filter 

propagation  cycle,  the  estimates  x,(t]  ,)  and  ,)  are 

I  c~  I  L~  J 

used  to  generate  control  signals  to  point  the  FLIR  sensor 
optica]  centerline  at  the  target. 

3.2.2  Measurement  Model.  As  an  alternative  to  the  64- 
dimensional,  non-linear  measurement  model  of  Equation  ( 2- 
41),  Rogers  [15]  developed  an  enhanced  correlation  algorithm 
to  provide  two  "measurements"  of  centroid  offsets  in  the 
FLIR  plane  to  a  linear  Kalman  filter  measurement  model.  The 
correlation  algorithm  enhancement  occurs  in  several  ways 
[18]. 

First,  the  current  FLIR  data  frame  is  correlated  with  a 
template  (an  estimate  of  the  target's  intensity  function), 
instead  of  being  correlated  with  the  previous  FLIR  data 
frame.  This  template  will  be  developed  in  Section  3. 2. 2.1. 
Second,  instead  of  outputting  the  peak  of  the  correlation 
function,  the  enhanced  correlator  outputs  the  center  of  mass 
of  that  portion  of  the  correlation  function  that  is  greater 
than  some  predetermined  lower  bound,  a  technique  known  as 
"thresholding".  Therefore,  the  enhanced  correlator  does  nor 
suffer  the  problem  of  distinguishing  global  peaks  from  local 
peaks,  as  do  many  conventional  "peak-finding"  correlation 
algorithms.  Third,  by  using  the  enhanced  correlation  al¬ 
gorithm,  the  FLIR/laser  pointing  commands  are  generated  via 
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the  Kalman  filter  propagation  cycle  rather  than  as  the  "raw 
measurement"  output  of  a  standard  correlation  algorithm. 
Finally,  the  Kalman  filter  state  estimate  is  used  to  center 
the  template,  so  that  the  offsets  seen  in  the  enhanced 
correlator  algorithm  should  be  smaller  than  in  the  conven¬ 
tional  correlator.  This  increases  the  amount  of  "overlap" 
between  the  actual  FLIR  data  and  the  stored  template,  thus 
improving  performance. 

The  outputs  of  the  enhanced  correlator  are  the  two  FLIR 
plane  x  and  y  offsets  of  Equations  (2-1)  and  (2-2).  It  is 
these  enhanced  or  "pseudo-measurements",  rather  than  the  raw 
FLIR  data,  which  are  used  in  the  linear  Kalman  filter  mea¬ 
surement  update  cycle.  An  overview  of  the  enhanced  corre¬ 
lation  process  is  presented  next.  For  a  more  detailed 
development,  see  Rogers  [15]. 

3. 2. 2.1  Template  Generation.  The  template  de¬ 
scribed  in  this  section  is  part  of  an  overall  data  process¬ 
ing  algorithm,  shown  in  Figure  3.1.  The  following  discus¬ 
sion  of  the  template  refers  to  the  "Form  Smoothed  Template" 
block  in  Figure  3.1. 

The  template  fan  estimate  of  the  target's  intensity 
function)  is  generated  by  averaging  the  N  most  recent  cen¬ 
tered  intensity  functions  observed  by  the  FLIR  sensor.  The 
intensity  functions  are  centered  on  the  FLIR  plane  by  use  of 
the  "shifting"  property  of  the  Fourier  transform,  since  the 
frequency  domain  is  where  the  correlation  is  being 
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determined.  The  value  of  the  memory  size  N  is  determined  by 
the  dynamics  of  the  target  intensity  function,  with  more 
dynamically  changing  functions  requiring  smaller  values  of 

N. 

For  online  filter  applications,  the  use  of  true  finite 
memory  filtering  techniques  causes  difficulties  with  con¬ 
straints  on  computer  memory.  To  avoid  these  difficulties, 
"exponential  smoothing"  was  used  to  approximate  the  averag¬ 
ing.  The  properties  of  exponential  smoothing  are  very 
similar  to  finite  memory  filtering  [7],  but  with  the  ad¬ 
vantage  of  only  requiring  storage  for  the  FLIR  data  from  the 
previous  sample  instead  of  N  previous  samples.  The  exponen¬ 
tial  smoothing  algorithm  maintains  the  template  by  the  use 
of  the  following  equation: 

=  vK*j)  +  (l-y)l(^;)  (3-11) 

where : 

A 

I (£  )  =  "smoothed  estimate"  (template)  of  target's 
intensity  function 

I(£-)  =  "raw"  intensity  function  from  the  current  FLIR 
data  frame 

y  =  smoothing  constant;  0  <  y  s  1 

The  smoothing  constant's  value  is  eomparaole  to  the  value 
selected  for  N.  Larger  smoothing  constant  values  emphasize 
the  current  FLIR  data  and  correspond  to  small  values  of  N. 
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Based  on  previous  studies  [3,5,11,13,14,17],  a  value  of  0.1 
will  be  used  for  the  smoothing  constant. 

Figure  3.1  details  the  structure  of  the  overall  enhanced 
correlation/linear  measurement  model  data  processing  al¬ 
gorithm.  This  structure  deals  with  an  algorithm  which  uses 
data  from  a  non-rotated  FLIR  field  of  view.  For  this  rhesis 
work,  the  algorithm  was  modified  to  use  a  diagonally  rotated 
FLIR  field  of  view.  That  modification  will  be  presented  in 
Chapter  4.  knen  the  incoming  FLIR  data  is  received,  a 
transformation  to  the  Fourier  domain  by  a  fast  Fourier 
transform  is  done.  The  FLIR  image  data  is  then  centered  in 
the  FLIR  plane  by  use  of  the  shifting  property  of  the  Four¬ 
ier  transform.  The  filter  computed  shift  is  equal  to: 

)  (  J-±2  ) 

K  =  vt{t)  +  ya(t()  (3-13) 

where  ,  xa,  ,  and  v ^  are  the  estimates  of  the  states 
defined  in  Equation  (3-1). 

There  are  several  reasons  for  performing  a  Fourier 
transformation  on  the  FLIR  data.  First  is  thai  the  correla¬ 
tion  described  in  Section  3. 2. 2. 2  is  readily  done  in  the 
Fourier  domain.  The  second  is  that  this  form  of  the  al¬ 
gorithm  allows  for  implementation  of  optical  processing 
i  15]  . 


* 
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The  shifts  of  Equations  (3-12)  and  (3-13)  (the  "X-i 
Phase  Shift"  block  of  Figure  3.1)  are  performed  using  the 
shifting  property  of  Fourier  transforms.  This  property 
allows  translation  shifts  in  the  spatial  domain  to  be  accom¬ 
plished  by  performing  a  linear  phase  shift  in  the  Fourier 
domain  [15].  The  phase  shift  is  expressed  as: 

/'<  g(  f.  ,  y-y^_ ,  J  I  =  G(/A;,/y)exp<  -j2x  l  . )}  (3-14) 

where : 

q(x,y)  -  2-dimensional  spatial  data  array 
/,  ■}  =  Fourier  transform  operator 
G(fx,fy)  =  F(q(x,y)  j 

After  centering,  the  data  is  incorporated  into  the 
template  using  the  exponential  smoothing  technique  of  Equa¬ 
tion  (3-11).  It  is  this  template  which  is  now  stored  and 
used  for  correlation  with  the  next  FLIR  data  frame  to  pro¬ 
duce  the  "pseudo-measurement"  used  by  the  filter. 

3 . 2 . 2 . 2  Enhanced  Correlation  "Pseudo-Measurements" . 
The  template,  as  developed  in  the  previous  section,  serves 
as  the  filter's  best  estimate  of  the  shape  of  the  target 
intensity  function  prior  to  receiving  a  new  FLIR  data  frame. 
The  correlation  of  the  incoming  FLIR  data  with  the  template 
provides  the  position  offset  of  the  target  intensity  func¬ 
tion  from  the  center  of  the  field  of  view.  This  cross-cor¬ 
relation,  performed  in  the  Fourier  domain  space,  is  computed 
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by  taking  the  inverse  fast  Fourier  transform  of  the  follow¬ 
ing  equation  [15]: 


/  :g(x,  V)  ■  l(x,y)  )  =  G(  t:< ,  )v )  L  (  tx,Sv)  (3-15) 

where : 

F{  i  =  Fourier  transform  operator 
q(  x ,  y  )  ■  l(x ,  y)  --  cross  correlation  of  q(x,y)  and  1  (x,y) 
g(x,y)  =  measured  target  intensity  function 
1  (x,y)  =  expected  target  intensity  function  (template) 
G(fx,fy)  =  F\z{x,y)  ( 

1>' {  f  K ,  f  v)  =  complex  conjugate  of  /{l(x,i/)> 

After  the  inverse  fast  Fourier  transformation  has  been  done, 
the  values  of  the  correlation  function  q{x , y)  ■  I(x ,  y)  are 
modified  so  that  any  value  in  the  correlation  function  less 
than  0.3  [14]  of  the  function's  maximum  value  is  set  to 
zero.  This  "thresholding"  technique  is  used  to  eliminate 
false  peaks  in  the  correlation  function  that  occur  due  to 
noise  and  other  effects.  As  indicated  in  Figure  3.1,  the 
output  of  the  inverse  fast  Fourier  transform  is  the  of i set 
of  the  "thresholded"  FLIR  intensity  function  from  the  center 
of  the  FLIR  field  of  view.  This  offset  is  assumed  to  be  the 
result  of  the  summed  effects  of  target  dynamics,  atmospheric 
jitter,  and  measurement  noise.  Expressed  in  terms  of  the 
filter  states  of  Equation  (3-1),  the  offset  measurement  in 
p i x  e 1 s  is: 
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These  two  measurements  can  then  be  represented  in  state 
space  form  as: 

z(^)  =  )  +  v/(</)  (3-18) 


where : 


z  (t.) 

x/(^) 

»/ 


f  xoffset '  y onset  - 

state  vector  of  Equation  ( 3-1  ) 

2x6  measurement  matrix 

2-dimensional,  discrete-t ime ,  zero-mean,  white 
Gaussian  measurement  noise  of  covariance  R; 


The  measurement  matrix,  H j,  is 


rioooio  "J 

H/  "  Lo  1  0  0  0  1  J 

and  the  covariance  matrix  Ry-  is  [9,15]: 

r  0.00363  0  -j 

R/  “  L  o  0.00598  J 


(  3-19) 


(  3-20  ) 


The  measurement  vector  and  measurement  distribution 
matrix  are  then  used  in  the  linear  Kalman  filter  update 
equations.  It  is  assumed  the  reader  already  has  a  knowledge 
of  linear  Kalman  filter  theory,  so  the  equations  will  not  be 
repeated  here. 


j . j  Center  of  Mass  Offset  Filter 

The  main  thrust  of  this  thesis  work  is  to  consider  the 
use  of  additional  measurements  to  aid  in  determination  of 
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the  missile  hardbody  location.  The  purpose  of  the  second 
filter  developed  in  this  effort  is  to  estimate  the  offset 
distance  between  the  filter  image  intensity  centroid  es¬ 
timate  and  the  missile  hardbody  center  of  mass.  The  deci¬ 
sion  [8]  was  made  to  orient  the  offset  distance  angularly 
using  the  filter  estimated  FLIR  plane  image  centroid  velo¬ 
city  vector.  This  decision  was  based  on  the  desire  to  use 
any  information  already  available  to  aid  in  determining  the 
missile  center  of  mass  location.  The  six-state  filter 
already  provides  estimates  of  the  FLIR  image  velocity,  and 
this  should  be  reflective  of  the  missile  hardbody  velocity 
vector.  Figure  3.2  shows  the  physical  representation  of  the 
filter's  estimate  of  the  offset  distance. 


3.3.1  Dynamics  Model.  For  this  effort,  the  offset 
between  the  filter  FLIR  image  centroid  and  the  missile 
center  of  mass  was  represented  as  a  bias.  The  decision  was 
made  [8]  to  represent  this  offset  bias  using  a  model  of  an 
integrator  driven  by  white  Gaussian  noise  (actually  pseudo¬ 
noise  for  filter  tuning:  an  undriven  integrator  yields  a 
bias  as  an  output).  The  single  state  representation  of  the 
linear,  time- invar iant ,  stochastic  differential  eguation  is 


where : 

x  ,  .  =  state  representinq  offset  distance  between 

missile  cm  and  FLIR  image  centroid 
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wcm  0fj.,P(  ( f- )  =  zero  mean,  white  Gaussian  noise  with 
strength  Q 

The  scalar  representation  of  the  filter  propagation 
Eguations  (3-6)  and  (3-7)  are: 


offset^  t-1  > 


>ff  set  ^  l  ^ 


(  3-22 


cn 7  offset 


c/t?  offset 


~d  cm  offset 


(  3-23 


3.3.2  Measurement  Model. 

One  of  the  reasons  for  the  method  in  which  the  center  of 
mass  offset  measurement  simulation  was  developed  in  Chapter 
2  was  to  preserve  linearity  in  the  filter  equations.  Other 
formulations  (such  as  separately  estimating  the  x  and  y 
components  of  the  bias)  were  considered,  but  involved  non¬ 
linear  measurement  models.  The  use  of  a  single  state,  repre¬ 
senting  the  magnitude  of  the  bias,  allowed  a  linear  filter 
formulation  to  be  used.  This  allows  the  scalar  measurement 
model  to  be  expressed  as: 


2  cm  offset^  xcm  offset  +  Vt  m  offset  ^  (3-. 

where : 

z  m  nff.:,,t  ( i )  =  the  niass  center  offset  measurement 

at  time  l1 

v  {(, )  =  discrete-time,  zero  mean,  white  Gaussian 

measurement  noise  with  variance  R 

3 . 4  Summary 

This  chapter  has  presented  moaels  upon  which  the  two 
Salman  filters  used  in  this  thesis  effort  are  based.  Ihe 
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six-state  filter  includes  models  to  estimate  target  dynamics 
and  atmospheric  jitter.  The  measurement  model  included  the 
development  of  "pseudo-measurements",  which  are  created  by 
correlating  the  current  FLIR  data  frame  with  a  template 
representing  the  target's  expected  infrared  intensity  func¬ 
tion  shape.  The  second  single  state  filter  provides  es¬ 
timates  of  the  offset  of  the  missile  hardbody  center  of  mass 
from  the  six-state  filter  estimated  FLIR  image  centroid. 

The  offset  was  modeled  as  bias,  the  output  of  an  undriven 
integrator,  and  the  white  driving  noise  was  added  to  allow 
for  filter  tuning.  The  offset  distance  is  oriented  using 
the  six-state  estimate  of  the  FLIR  plane  FLIR  image  velocity 
vector.  The  measurement  model  for  this  filter  presented  the 
measurement  as  a  direct  representation  of  the  offset  state 
plus  measurement  noise. 
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4 .  Tracking  Algorithm  and  Performance  Evaluation  Tools 

4 . 1  Introduction 

The  purpose  of  this  chapter  is  to  present  an  overall 
description  of  the  tracking  algorithm  used  in  this  thesis 
effort.  Modifications  to  the  material  presented  xn  Chapter 
3  to  accommodate  a  diagonally  rotated  FL1R  f ield-of -view 
will  be  presented.  The  truth  model  and  filter  model  para¬ 
meters  used  in  this  work  will  also  be  detailed.  Finally, 
the  evaluation  tools  for  measuring  filter  performance  will 
also  be  shown. 

4 . 2  Tracking  Algorithm  Overview 

The  main  objective  of  this  thesis  research  is  to  develop 
an  algorithm  which  can  accurately  track  the  center  of  mass 
of  a  ballistic  missile  (in  boost  phase)  using  both  passively 
(FLIR)  and  actively  (low  power  laser  reflection  off  missile 
hardbody)  acquired  measurements.  This  involves  development 
of  two  separate  filters  and  of  the  simulation  of  the  active¬ 
ly  acquired  measurement. 

Figure  4.1  details  the  overall  algorithm.  The  raw  FLIR 
measurement  is  input  into  the  enhanced  correlator  data- 
processing  algorithm.  This  algorithm  is  contained  within 
the  dotted  lines  of  Figure  3.1  in  the  previous  chapter. 

After  the  "pseudo-measurements"  are  produced,  they  become 
linear  measurements  entered  into  the  six-state  linear  Kalman 
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filter.  The  filter  produces  estimates  of  the  FLIR  plane 
FLIR  image  centroid  position  and  velocity  vector.  These 
estimates  are  used  by  the  mass  enter  measurement  algorithm 
to  produce  a  "pseudo-measurement”  representing  the  output  of 
an  optical  sensor  receiving  reflected  laser  light  from  the 
missile  hardbody  as  the  laser  is  swept  up  and  down  the 
velocity  vector  direction  from  the  estimated  centroid  loca¬ 
tion.  This  output,  in  the  form  of  an  offset  distance,  is 
then  used  by  a  second  filter  to  produce  an  estimate  of  the 
target  center  of  mass  position.  This  estimate  would  then  be 
the  primary  means  of  aiming  a  high  energy  laser  at  the 
missile  hardbody. 

4 • 3  Field-of-View  Rotation 

In  f ield-of-view  rotation  schemes,  the  FLIR  sensor  field 
of  view  is  rotated  so  that  the  target's  filter-estimated 
FLIR  plane  velocity  vector  is  aligned  with  one  of  the  FLIR 
plane  axes  or  with  some  other  appropriate  direction  in  the 
FLIR  plane.  The  ides,  first  investigated  by  Leeney  [4;,  was 
an  attempt  to  maintain  lock  on  highly  dynamic  targets  that 
could  "jink"  in  either  of  the  two  FLIR  plane  directions. 
Norton  [13]  extended  Leeney's  preliminary  work  by  examining 
the  use  of  a  rotating  f ield-of-view  in  an  MMAF  scheme.  Most 
recently,  Rizzo  [14]  investigated  the  performance  of  identi¬ 
cal  filters  with  different  FOV  rotation  schemes.  He  com¬ 
pared  non-rotating,  rotating  (aligning  the  principal  axis  of 


the  f ield-of-view  with  the  estimated  velocity  vector),  and 
diagonally  rotating  (aligning  the  diagonal  of  the  field-of- 
view  with  that  estimated  velocity  vector)  fields  of  view. 
The  filter  using  the  diagonally  rotated  f ield-of-view 
( DRFOV )  had  the  lowest  mean  tracking  error  and  standard 
deviation.  Based  on  Rizzo's  [14]  work,  it  was  decided  to 
use  a  DRFOV  in  conjunction  with  the  six-state  filter. 

As  previously  mentioned  in  Chapter  3,  the  use  of  a 
rotating  field  of  view  causes  some  modifications  to  be  made 
to  the  data  processing  algorithm  of  Figure  3.1  of  the  pre¬ 
vious  chapter.  These  modifications  are  shown  in  Figure  4.2 
with  the  addition  of  the  "rotate"  blocks.  These  modifica¬ 
tions  will  be  described  briefly  here;  for  a  complete  devel¬ 
opment  see  Rizzo  [14]. 

The  basis  for  the  rotating  FOV  is  the  filter  estimated 
target  FLIR  plane  velocity  vector  orientation.  This  is 
given  by: 


(  4-1  ) 

where  Qr  is  shown  in  Figure  3.2 

Note  that  the  terms  in  Equation  (4—1 )  are  the  third  and 
fourth  states  in  the  six-state  filter  defined  in  Equation 
i  3-1).  This  allows  the  filter  to  provide  the  orientation 
control  directly  to  the  FLIR  sensor  for  on-line  application. 
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The  negative  sign  in  the  numerator  of  Eguation  (4-1)  keeps 
the  orientation  angle  positive  in  the  counterclockwise 
direction  from  the  FLIR  plane  x  axis  when  viewed  from  the 
inertial  frame  origin.  This  also  allows  a  direct  comparison 
with  the  truth  model  velocity  orientation  angle,  facilitat¬ 
ing  the  calculation  of  the  error  statistics. 

To  simulate  the  physical  rotation  of  the  FLIR  sensor, 
the  incoming  FLIR  data  is  rotated  before  entering  the  data 
processing  algorithm  of  Figure  4.2.  This  is  simulated  by 
performing  a  negative  rotation,  based  on  the  positive  orien¬ 
tation  angle,  on  the  location  and  orientation  of  the  in¬ 
dividual  Gaussian  intensity  functions  described  in  Section 
2.4.1.  This  corresponds  to  a  positive  rotation  of  the  FOV 
and  aligns  the  FOV  with  the  positive  orientation  angle. 
Mathematically,  this  is  done  by  first  rotating  the  intensity 
function  peaks  defined  by  Eguation  (2-40)  using  the  follow¬ 
ing  rotational  transformation  matrix,  first  produced  by 
Norton  [!_}]: 


“  y  '  “ 

peak 

~  A 

COS0  J 

A 

-sin Qj 

x  peak 

_  Y  peak  _ 

A 

sin 

A 

COS0  j 

V 

■  peak  _ 

Note  here  that  the  primed  variables  are  in  the  rotated 
coordinate  system.  The  intensity  function  of  Eguation  ( 2- 
40)  (refer  to  Figure  2.4)  then  becomes: 


11*'  'J/'  •*' oeak^'  ’V'wkV)  i  =  ^  maxexP  i  -0 . 5  [  Aa.  '  Ay  '  : 


-'ea/c 


P  '  A.'VAi/'  u  \ 


<  4  —  3  ) 
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where : 


A*'  =  (X*  -x' peafr )  cosA0  +  (  y’  ~y' peak  )  sinA© 

Ay'  =  0//-y';)eaA)cosA©  -  (*'-*' peak)  sinA9 

A©  =  difference  between  the  truth  model  velocity 
orientation  angle  and  the  filter  computed 

A 

orientation  angle,  i.e.  A©  =  Qj  -  Q j 
y' t  y'  =  rotated  coordinates  from  the  original  FLIR 
coordinate  frame  via  Equation  (4-2) 

Once  the  incoming  FLIR  data  has  been  rotated,  the  data 
processing  algorithm  generates  the  template  in  the  same 
manner  as  was  done  in  Figure  3.1. 

Recall  that  the  incoming  data  was  centered  in  the  FLIR 
plane  via  the  shifting  of  Equations  (3-12)  and  (3-13). 
However,  these  shifts  were  computed  in  the  unrotated  filter 
coordinate  system,  while  the  current  FLIR  data  image  is  in  a 
rotated  frame.  To  implement  the  algorithm  properly,  the 
translational  shift  is  accomplished  in  the  rotated  frame  by 
the  transformation: 

*’ shift  1  =  cos®/  -sin©/  r  x shin 

i/'  .  a  a  I  y  .  .. 

_  v  shift  J  sm©^  cosQf-  I  L  shift 

where  x^hift  anc*  Vrhift  are  9^ven  Equations  (  3 
This  transformation  is  accomplished  by  the 


block  following  the  tracking  algorithm  block. 


With  rhis 


transformation,  both  the  current  image  data  and  the  template 
are  in  the  same  coordinate  frame. 


The  final  modification  to  the  algorithm  occurs  with  the 
"rotate"  block  after  the  " IFFT"  block.  The  outputs  of  the 
IFFT  block  are  the  linear  offsets  between  the  current  image 
data  and  the  centered  template.  However,  these  offsets  are 
in  a  rotated  frame,  while  the  states  in  the  filter  are  in 
the  original  unrotated  frame.  To  insure  compatibility,  the 
offsets  are  rotated  back  into  the  original  frame  by  the 
transformation : 
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(4-5) 


where  the  primed  and  unprirced  coordinates  are  in  the 
unrotated  and  rotated  coordinate  systems,  respectively. 


4 . 4  Truth  Model  Parameters 

Since  the  scenario  used  for  this  thesis  was  a  continua¬ 
tion  of  that  used  by  Rizzo,  the  truth  model  parameters  are 
the  same.  The  rationale  for  the  choice  of  these  parameters 
is  fully  detailed  by  Rizzo  [14];  it  will  not  be  repeated 
here.  For  the  nominal  ballistic  missile  trajectory  studied, 
the  initial  conditions  for  the  target  inertial  frame  (see 
Section  2.3.1)  position  and  velocity  vectors  are: 

ev  -  2Q,noO  meters 
e,,  =  10  0,000  meters 
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ev  =  2,000,000  meters 


v* 


-2500 


meters 

second 


4330 


meters 

second 


v,  = 


meters 

second 


The  maximum  intensity  value  for  each  of  the  intensity 
functions  of  Equation  (2-40)  is  20  intensity  units.  The  RMS 
value  of  v ,  which  is  the  sum  of  spatially  correlated 
background  noise  (t>jk)  and  FLIR  sensor  noise  (njk)  in 
Equation  (2-41),  is  one.  This  yields  a  signal-to-noise  ratio 
of  twenty,  which  is  typical  of  many  tracking  scenarios  [18]. 

In  his  investigation  of  plume  "pogo"  effects,  Rizzo 
[14]  desired  to  have  the  exhaust  plume  in  the  FLIR  plane 
contained  within  a  5  x  5  pixel  window.  This  value  was 
chosen  so  that  the  maximum  plume  "pogo"  oscillation  would 
"fit"  within  an  8  x  8  pixel  FOV.  To  represent  the  reference 
ellipsoid  hotspot  dispersion  of  Figure  2.4,  the  hotspot 
dispersion  along  the  eFV  direction  of  the  target  frame  (Equa¬ 
tion  2-35)  was  chosen  to  be  1  pixel,  and  the  dispersion  in 
the  e  direction  as  1.5  pixels.  The  pixel  proportionality 
constant  was  initially  continued  from  Rizzo's  value  of  1.5 
micro-radians/pixel.  Although  the  "pogo”  phenomenon  was  not 
considered  in  this  research,  the  desire  to  keep  the  missile 
hardbody  center  of  mass  and  the  missile  plume  centroid  in 


80 


the  same  8x8  pixel  FOV  prompted  the  continued  use  of  the 
same  values  to  describe  the  hotspot  dispersion. 

The  variance  and  mean  squared  value  for  the  atmospheric 
jitter,  in  both  x  and  y  directions,  was  continued  from 
previous  theses  [13,14,18]  as  0.2  pixels*'.  The  truth  model 
bending/vibration  parameters  are  in  Section  2.2.3. 

4 . 5  Filter  Parameters 

For  the  first  analysis  of  the  problem  involving  use  of 
actively  acquired  measurements,  it  was  decided  to  concen¬ 
trate  solely  on  the  problem  of  tracking,  rather  than  acqui¬ 
sition  and  tracking.  Therefore,  in  the  simulation  used  for 
this  research,  the  six-state  filter  is  initialized  with  zero 
errors  in  the  position  and  velocity  at  time  t=0.  The  atmo¬ 
spheric  jitter  states  are  also  initialized  to  zero.  The 
initial  state  covariance  matrix  is  identical  to  that  used  by 
Rizzo,  for  states  which  are  common  to  both  filters.  with 
this  condition,  the  matrix  becomes: 


P(t  ) 


“  10  00  0  0  0“ 
0  10  0  0  0  0 

0  0  2000  000 
000  2000  0  0 
0  0  0  0  0.2  0 

_000  0  0  0 . 2  _ 


(4-6) 
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For  the  initial  conditions  on  the  single  state  filter, 
both  the  state  and  the  covariance  were  set  to  zero.  It  must 
be  noted  that  using  these  initial  values  will  have  a  signif¬ 
icant  effect  on  the  initial  transient  performance  of  the 
filter.  The  values  for  the  dynamics  and  measurement  noise 
strength  were  the  subjects  of  analysis  in  both  filters,  and 
these  will  be  discussed  in  Chapter  5. 


4 . 6  Tracking  Algorithm  Statistics 

The  performance  of  the  tracking  algorithm  is  evaluated 
using  Monte-Carlo  simulation  techniques  [6].  Previous 
research  has  shown  that  ten  Monte-Carlo  runs  demonstrate 
sufficient  convergence  to  the  actual  statistics  obtained 
fron  an  infinite  number  of  runs  [17,18].  Based  on  these 
previous  efforts,  ten  Monte-Carlo  runs  were  used  to  analyse 
tracker  performance  in  this  research. 

The  sample  mean  errors  of  the  tracking  algorithm's 
estimates  are  calculated  as  [17]: 


«,■>  - 


|  I[s 

«=: 


(*/)  - 


Tn'U 


where : 


(4-7) 


F.  s{t  ■)  =  sample  mean  error  of  state  estimate  at  time  £  ■ 
averaged  over  V  runs 

W* t>  =  state  estimate  at  time  t,  during  simulation  n 
sTn^ti')  =  truth  model  state  value  at  time  t{  during 
simulation  n 

V  =  number  of  Monte-Carlo  runs 
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The  sample  variance  of  the  error  is  given  by: 


a 


s 


_1_ 

•V- 1 


n= . 


In 


V- 1 


(4-8) 


These  two  statistics  are  calculated  for  the  target  x  and 
y  positions  and  velocities  due  to  dynamics,  and  the  FLIR 
image  centroid  x  and  y  coordinates.  The  performance  of  the 
filter  in  estimating  the  centroid  location,  as  well  as 
velocity  is  of  primary  importance  ’ n  generating  the  mass 
center  offset  measurement  used  by  the  single  state  filter. 
The  same  error  statistics  are  also  calculated  for  the  offset 
estimate.  All  the  error  statistics  are  in  units  of  pixels, 
and  calculated  for  both  before  and  after  measurement  update. 

For  further  simplification,  the  above  statistics  are 
temporally  averaged  over  the  seven  second  simulation.  These 
temporal  averages,  when  used  in  conjunction  with  the  plotted 
time  histories  of  the  error  statistics,  provide  some  indica¬ 
tion  of  trends  during  the  simulation. 


4 . 7  Performance  Plots 

Thirteen  plots  are  used  to  assess  filter  performance  in 
this  study.  The  first  ten  plots  are  used  to  describe  the 
performance  of  the  six-state  f  "  r.  They  are: 

1.  True  *  position  rms  error  vs.  filter-computed  x  position 
rms  error 

2.  True  v  position  rms  error  vs.  filter-computed  y  position 
rms  error 
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.  .  4- 

3.  Mean  *  target  position  error  ±  one  a  at  all  time  l- 

4.  Mean  v  target  position  error  ±  one  a  at  all  time  t- 

5.  Mean  a;  target  position  error  ±  one  a  at  all  time  t  ■' 

6.  Mean  y  target  position  error  ±  one  a  at  all  time  t- 

7.  Mean  *  centroid  position  error  ±  one  a  at  all  time  t- 

8.  Mean  y  centroid  position  error  ±  one  o  at  all  time  t- 

9.  Mean  a:  centroid  position  error  ±  one  a  at  all  time  t~ 

10.  Mean  y  centroid  position  error  ±  one  a  at  all  time  t- 

Examples  of  plots  1,  3,  7,  9,  11,  and  13  are  shown  in 
Figures  4 . j  through  4.8.  The  first  two  plots  indicate  the 
adequacy  of  the  filter  tuning  by  direct  comparison  of  the 
actual  true  rms  error  vs.  the  filter  computed  rms  error. 

Note  that  in  Figure  4.3,  the  constant  graph  is  the  filter 
computed  error  and  the  graph  with  the  peaks  and  valleys 
represents  the  actual  error.  The  degree  of  overlap  between 
the  two  graphs  indicates  how  well  the  filter  is  tuned 
against  the  truth  model.  Plots  3  through  6  indicate  how  the 
tuning  is  affecting  the  filter's  ability  to  estimate  the 
portion  of  the  image  position  due  strictly  to  target  dyna¬ 
mics.  Plots  7  through  10  provide  primary  tracking  perfor¬ 
mance  information,  since  the  location  of  the  image  centroid 
is  used  in  the  determination  of  the  hardbody  center  of  mass. 
Plots  11-13  assess  the  performance  of  the  filter  which 
estimates  the  offset  distance  between  the  FLIR  image 
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T  X  ,  T  Y 


Figure  4.3  Example  of  Actual  vs.  Filter  Computed  RMS  Error 
Performance  Plot  (Plot  #1) 
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IMK  IN  SECOND 


T  X  ,  T  Y 


U1  ® 

a 


M  a;  a;  o  a:  —  z:  x  —  x  M  —  zr. 


Figure  4.4  Example  of  Position  Error  at  time  t 


Performance  Plot  (Plot  #3) 


ESTIMATED  X  PLUS  CENTROID  POS 


Performance  Plot  (Plot  #9) 


centroid  and  the  missile  hardbody  center  of  mass.  The  error 
for  this  filter  is  determined  in  two  steps.  First,  the 
filter  er -imated  location  of  the  center  of  mass  is  compute! 
using  the  filter  estimated  offset  length,  the  estimated 
velocity  vector  orientation  angle,  and  the  estimated  image 
centroid  location.  The  error  is  then  calculated  as  the 
difference  between  the  estimated  center  of  mass  and  the 
truth  model  center  of  mass.  As  previously  mentioned  in 
Section  2.4.2,  there  may  be  sample  times  in  the  simulation 
when  no  measurement  occurs.  In  this  case,  the  last  filter 
estimate  based  on  a  measurement  update  is  used  for  error 
computation  until  the  next  measurement  update  occurs.  The 
filter  error  plots  are: 

11.  True  offset  rms  error  vs  filter  offset  rms  error 

12.  Mean  offset  error,  ±  one  a,  at  all  time  #■ 

13.  Mean  offset  error,  ±  one  a,  at  all  time  t~ 

Similar  to  plots  1  and  2,  plot  11  indicates  the  adeguacy  of 
filter  tuning.  Plots  11  and  12  are  the  primary  indicators 
of  filter  performance,  indicating  the  filter's  ability  to 
locate  the  missile  center  of  mass. 

4.7.1  Plot  Designation  Codes.  The  first  ten  plots  have 
four  fields  in  the  graph  label.  The  first  field,  labeled 
"TX,  TY" ,  indicates  the  correlation  time  constant  (in 
seconds)  used  for  the  velocity  model  in  the  x  and  y 
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directions.  The  second  field,  labeled  "TAF" ,  is  the  time 
constant  (in  seconds)  usea  for  the  atmospheric  jitter  model, 
and  is  equal  in  both  x  and  y  directions.  The  third  field, 
labeled  "VD" ,  is  the  variance  of  the  discrete  driving  noise 
for  the  velocity  states  in  pixels2/seconds‘,  and  is  assumed 
equal  in  both  the  x  and  y  directions.  The  fourth  field, 

"VA" ,  is  the  variance  of  the  driving  noise  associated  with 
the  atmospheric  jitter,  and  is  also  assumed  equ-i  in  both 
the  x  and  y  directions.  It  also  has  units  of  pixels  - 
/seconds2.  The  last  field,  "BND",  simply  indicates  the 
inclusion  of  bending  states  in  the  truth  model. 

For  the  three  plots  associated  with  the  offset  filter, 
there  are  five  fields  in  the  graph  label.  The  first,  "FLT 
DN" ,  indicates  the  variance  of  tne  filter  discrete-time 
dynamics  driving  noise,  in  pixels2.  The  next  field,  "MN" , 
indicates  the  variance  of  the  filter  measurement  noise,  also 
in  pixels2.  The  third  field,  "  TR  MN" ,  indicates  the  mea¬ 
surement  noise  variance  the  truth  model  generation  of  the 
offset  measurement  used  by  the  filter.  The  fourth  field, 
"OFF",  is  the  true  offset  distance,  in  meters,  between  the 
missile  center  of  mass  and  the  truth  model  image  centroid. 
The  last  field,  "PC",  indicates  the  pixel  proportionality 
constant,  in  micro-radians/pixel.  The  values  given  to  these 
parameters  will  be  discussed  in  Chapter  5. 
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4 . 8  Summar 


This  chapter  has  presented  the  overall  tracking  algo¬ 
rithm  used  in  this  research  by  combining  the  results  of 
Chapters  2  and  3.  An  overall  view  of  the  algorithm  was 
presented,  along  with  modifications  to  the  algorithm  caused 
by  the  inclusion  of  a  diagonally  rotated  FLIR  field  of  view. 
The  truth  model  and  filter  parameters  used  in  this  research 
were  then  detailed.  Finally,  the  statistical  tools  used  for 
the  filter  performance  evaluation  were  shown. 
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5. 


5 . 1  Introduction 

This  chapter  presents  the  results  of  the  evaluation  of 
the  performance  character istics  of  the  algorithms  discussed 
in  Chapter  4.  Section  5.2  details  the  observability  analys¬ 
is  conducted  to  determine  the  states  which  were  used  in  the 
six-state  filter.  Section  5.3  then  presents  the  analysis  of 
the  tuning  process  for  that  filter.  Once  tuning  of  the  six- 
state  filter  was  accomplished,  a  parameter  study  involving 
the  single  state  filter  for  estimating  the  mass-center-to- 
intensity-centroid  offset  was  undertaken.  Section  5.4 
presents  the  results  of  the  sensitivity  study,  in  which  the 
parameters  of  truth  model  and  filter  measurement  noise 
variance  were  altered.  Section  5.5  shows  the  corresponding 
robustness  study  results,  in  which  the  filter  was  not  in¬ 
formed  of  parameter  changes  in  the  real  world.  Section  5.6 
and  5.7  then  pr  sent  the  results  from  the  parameter  studies 
involving  the  truth  model  offset  distance  and  pixel  propor¬ 
tionality  constant. 

5 . 2  Filter  Observability  Study 

Rizzo  [14],  in  his  research,  discovered  an  observa¬ 
bility-like  problem  in  the  filter  models  used  to  describe 
the  "pogo"  phenomenon  he  was  investigating.  He  compared  the 
performance  of  an  eight-state  benchmark  filter  (modelling 
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FLIR  plane  image  position,  velocity,  acceleration,  and 
atmospheric  jitter  in  both  x  and  y  directions)  to  a  ten- 
state  filter  (modelling  the  same  eight  states  plus  the  two 
"pogo"  states,  pogo  effect  position  displacement  and  velo¬ 
city)  which  also  included  modelling  for  "pogo"  effects.  The 
performance  of  both  filters  was  measured  against  a  14-state 
truth  model  which  also  included  "pogo"  effects  modelling 
(the  same  12-state  model  described  in  Section  2.2  plus  the 
two  "pogo"  states).  In  examining  the  results,  he  discovered 
a  performance  anomaly,  in  that  the  eight-state  filter  out¬ 
performed  the  ten-state  filter.  Since  the  higher  order 
filter  had  "knowledge"  of  the  "pogo"  effects  in  the  truth 
model  and  the  lower  dimensioned  filter  did  not,  this  was  the 
opposite  of  the  results  to  be  expected.  In  investigating 
this  anomaly,  he  performed  a  stochastic  observability  test 
[6]  on  both  the  eight-  and  ten-state  filters.  The  velocity 
and  acceleration  states  were  found  to  be  weakly  observable, 
particularly  the  acceleration  states.  Rizzo  [14:7-5]  recom¬ 
mended  that  alternate  acceleration  models  be  considered  in 
the  future,  or  that  the  acceleration  states  be  dropped  from 
the  filter  completely,  in  order  to  address  this  difficulty. 

For  this  research,  it  was  decided  [8]  to  reduce  the 
filter  order  by  dropping  acceleration,  and  possibly  velocity 
states,  from  the  eight-state  filter.  To  determine  the  final 
order  of  the  filter,  the  observability  was  reexamined.  The 
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stochastic  observability  condition  to  be  examined  is  given 


by  the  following  relationship  [6:243]: 

<*I  ^  I  '-(t  -mt. )  <  81  (5-1) 

j~  i  ;V~  /  J  1  1  1  1 

where  the  summation  term  in  Equation  (5-1)  is  the 
stochastic  observability  Grammian  matrix.  If  there  exist 
positive  numbers  a  and  6,  and  0  <  a  <  B  <  oo  f  and  a  posi¬ 
tive  integer  ,V  such  that,  for  all  i  >  N,  the  above  relation¬ 
ship  holds,  then  the  system  is  said  to  be  stochastically 
observable.  Because  of  numerical  precision  problems 
involved  with  exponential  terms  within  the  Grammian  matrix, 

.V  =  13  and  i  =  14  were  used  in  the  test. 

The  observability  Grammian  matrix  (diagonal  terms)  for 
Rizzo's  eight-state  benchmark  filter  is: 

“  3856  -  --  -- 

-  2341  -  -  -  - 

-  310  -  -  - 

-  -  -  188  - 

-  -  -  11.2  - 
-----6.8 

-----  244146644 

_  _____  -  148202729  _ 

The  order  of  the  states  in  the  observability  Grammian 
matrix  is:  position  (x  and  y),  velocity  (x  and  y),  accelera¬ 
tion  (x  and  y),  and  atmospheric  jitter  (x  and  y) .  It  is 
necessary  to  note  that  the  observability  matrix  was  not  a 
diagonal  matrix:  for  the  purpose  of  clarity  and  emphasis, 
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only  the  diagonal  terms  are  shown.  By  observation  of  the 
size  of  the  diagonal  terms,  states  with  potential  obser¬ 
vability  problems  can  be  distinguished  by  their  small  mag¬ 
nitude  relative  to  the  other  diagonal  entries.  For  determi¬ 
nation  of  system  model  observability,  the  Grammian  matrix 
eigenvalues  must  be  examined.  The  eigenvalues  for  the 
eight-state  filter  are: 

1673 
2756 
14.97 
24 . 66 
.03037 
.  0501 
148203577 
_ 244148042  _ 

Of  the  eight  eigenvalues,  the  third  through  the  sixth 
are  smaller  by  three  to  five  magnitudes  than  the  next 
smallest  eigenvalues.  Here,  eigenvalues  three  and  four 
correspond  basically  to  the  velocity  states,  and  five  and 
six  to  the  acceleration  states.  Particularly  the  accelera¬ 
tion  state  eigenvalues,  being  almost  zero,  indicate  almost 
complete  unobservability.  This,  in  combination  with  the 
relatively  small  diagonal  values  from  the  observability 
matrix,  led  to  the  decision  to  delete  the  acceleration 
states  from  the  model.  Even  though  the  velocity  states  were 
much  less  observable  than  either  the  position  or  atmospheric 
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jitter  states,  the  decision  was  made  to  keep  these  states  in 
the  model.  The  need  for  a  filter  estimate  of  the  velocity 
vector  orientation  angle  as  a  means  of  aiming  the  low  power 
laser  necessitated  this  decision.  With  the  deletion  of 
acceleration  from  the  eight-state  model,  and  the  remodelling 
of  velocity,  the  six-state  model  detailed  in  Chapter  3 
became  the  basis  for  the  filter  used  in  this  research  to 
estimate  FLIR  plane  target  parameters. 

5 . 3  Filter  Tuning  Process 

The  purpose  of  this  section  is  to  describe  the  tuning 
procedure  and  results  obtained  during  the  tuning  of  the  six- 
state  filter  for  the  benign  trajectory  used  for  this  re¬ 
search.  The  tuning  parameters  used  in  this  study  are  the 
variances  (in  both  the  x  and  y  directions)  of  the  first- 
order  Gauss-Markov  representations  of  the  target  velocity, 
and  the  variances  (in  the  x  and  y  direction)  of  the  atmo¬ 
spheric  jitter  position  processes.  This  decision  was  made 
after  conducting  an  initial  tuning  study  involving  the 
variances  and  correlation  times  of  both  the  target  velocity 
and  jitter  position  processes.  This  study  revealed  the 
sensitivity  of  the  filter  mainly  to  changes  in  the  varian¬ 
ces.  Additionally,  in  previous  research  [13,14,18],  a  value 
of  0.707  seconds  had  been  shown  to  be  representative  for  the 
atmospheric  jitter  correlation  time  constant,  and  the  deci¬ 
sion  [8]  was  made  to  remain  with  this  value  throughout  the 
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tuning  process.  A  value  of  8.5  seconds  was  used  [8]  for  the 
correlation  time  constant  of  the  velocity  process,  and  was 
also  kept  constant  throughout  the  tuning  process.  This 
value  was  selected  based  on  the  fact  that  the  dynamic  maneu¬ 
vering  expected  from  a  ballistic  missile  in  boost  phase 
would  be  somewhat  less  than  that  expected  from  a  large 
manned  bomber  (with  a  correlation  time  of  approximately  six 
seconds ) . 

Table  5.1  lists  the  variation  of  the  two  tuning  para¬ 
meters,  the  target  dynamics  velocity  and  jitter  position 
variances,  for  the  tuning  process.  For  the  first  tuning 
run,  the  dynamics  variance  was  set  at  5,  and  the  jitter 
variance  was  set  at  0.2.  The  time  averaged  statistics  are 
listed  in  Table  5.2  and  the  performance  plots  are  given  in 
Appendix  B,  Figures  B.1-B.10.  The  plots  in  Figure  B.l  and 
B.2  indicate  the  filter  is  underestimating  the  error  varian¬ 
ces  on  target  position  estimates  in  both  the  x  and  y  direc¬ 
tions.  The  dynamics  state  errors  in  Figures  B.2-B.6  also 
indicate  the  need  for  further  tuning.  The  centroid  plots, 
Figures  B.7-B.10,  are  very  important.  In  this  research,  the 
need  for  the  filter  to  provide  an  accurate  estimate  of  the 
FLIR  image  centroid  is  critical  to  the  aiming  of  the  low 
power  laser  and  resulting  estimation  of  the  missile  center 
of  mass  (  velocity  estimation  accuracy  is  also  important,  in 
order  to  establish  the  estimated  orientation  angle  of  Equa¬ 
tion  (4-1)  that  indicates  the  angular  orientation  of  the 
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center  of  mass  from  the  intensity  centroid).  Here,  par¬ 
ticularly  in  the  plots  for  the  centroid  position  error  after 
update  (Figures  B.9-B.10),  there  is  an  increasing  error  over 
time . 


Table  5.1.  Tuning  Run  Parameters 


RUN  # 

[  pixels' 

°x~>°v\ 

y  1  seconds 

j  a a2  ( pixels" 

1 

5 

0 . 2 

2 

50 

C  .  2 

3 

200 

0 . 2 

4 

800 

0 . 8 

Table  5.2.  Run  1  Temporally  Averaged  Statistics 


Error  in  Estimate  of : 

Mean 

1  Sigma 

x(t-) 

-1.4587 

.94001 

y{tt) 

x(t-) 

y(t-) 

2 . 5274 

1 .02790 

-1 . 2958 

. 90143 

2 . 2610 

.98574 

*c(V) 

-  .  0126 

.98152 

KC(V) 

1 .4844 

.93751 

3C(V> 

-  .  1263 

.  28028 

.  5154 

.  19054 

For  the  second  run,  the  target  dynamics  model  velocity 
variance  was  increased  to  50,  while  the  jitter  variance  was 
kept  at  0.2.  The  tuning  results  are  presented  in  Figures 
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Table  5.3.  Run  2  Temporally  Averaged  Statistics 


Error  in  Estimate  of : 

Mean 

1  Sigma 

*(ti  ) 

- . 56178 

1 . 05820 

v{tt) 

£(V) 

yit[) 

.75109 

1 . 10460 

-  .47145 

.  99714 

.61261 

1 .04780 

M*<  ) 

-.33741 

1 .00710 

vv  > 

.44683 

. 95115 

mV) 

-.  12905 

. 29073 

MV) 

.11751 

. 18595 

B.11-B.20  and  Table  5.3.  Although  some  improvement  in 
matching  the  filter-computed  error  variance  to  the  actual 
error  variance  has  been  made,  the  improvement  is  not  sig¬ 
nificant.  A  significant  improvement  in  the  target  dynamics 
position  state  error  and  especially  in  the  centroid  error 
has  been  achieved,  as  indicated  by  Figures  B.13-B.20.  In 
comparing  Figures  B.9  and  B.10  to  B.19  and  B.20,  note  the 
mean  error  line  slope  has  been  significantly  reduced  in  the 
x  channel.  The  y  channel  has  also  improved,  but  not  to  the 
same  degree.  This  is  due  to  the  fact  that  the  R  matrix  of 
Equation  (3-20)  in  the  measurement  model  weights  the  y  axis 
FLIR  plane  measurements  more  heavily  than  the  x  axis,  and 
thus  the  filter  y  axis  estimates  respond  less  to  changes  in 
the  dynamics  driving  noise  strength.  The  increase  in  the 
error  over  time  indicates  that  filter  divergence  would 
eventually  become  large  enough  to  cause  loss  of  tracking. 
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However,  a  target  reacquisition  algorithm  [18]  has  been 
previously  developed  in  past  work  as  part  of  the  simulation. 
If  the  errors  did  become  large  enough,  the  simulation  uses 
this  algorithm  in  a  target  acquisition,  rather  than  a  target 
tracking  mode.  This  would  prevent  total  divergence  from 
occurring . 

Several  interim  runs  were  performed,  in  which  the  target 
dynamics  velocity  strength  was  systematically  increased 
while  keeping  the  jitter  variance  constant.  A  final  value 
of  200  was  used  for  the  target  dynamics  velocity  variance. 
The  results  of  that  run  are  shown  in  Figures  B.21-B.30  and 
Table  5.4.  As  can  be  seen  from  Figures  B.21  and  B.22, 
changing  the  target  dynamics  velocity  variance  does  not 
significantly  affect  the  match  between  the  filter-computed 
error  variance  and  the  actual  error  variance.  However,  the 
errors  in  the  dynamics  states  and  particularly  in  the  cen¬ 
troid  location  have  been  significantly  reduced.  Note 
Figures  B.29  and  B.30.  The  error  in  the  x-axis  centroid 
location  is  essentially  zero.  The  y-axis  centroid  error, 
while  not  zero,  is  less  than  0.2  pixels  at  the  end  or  the 
simulation . 

The  next  step  was  to  discontinue  the  tuning  via  the 
target  dynamics  error  variance,  and  to  continue  the  tuning 
using  the  jitter  variance.  When  tuning  via  the  jitter 
variance,  a  counterproductive  result  was  noted.  As  the 
jitter  variance  was  increased,  the  match  between  the  filter- 
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Table  5.4.  Run  3  Temporally  Averaged  Statistics 


Error  in  Estimate  of: 

Mean 

1  Sigma 

A  , 

*(*t  ) 

- .  30019 

1 . 2293 

y(tt  ) 

A  ,  - 

. 53239 

1.1514 

*(*i  ) 

A  ,  - 

-  .23704 

1  .  1340 

y{tt  ) 

A  ,  ,  - 

.43096 

1 .0550 

) 

A  ,  , 

-  .  09359 

1 .0376 

vc^i  ] 

.  30553 

.98254 

*c(V) 

- .01949 

.31303 

. 12094 

. 18187 

computed  error  variance  and  the  actual  error  variance  in¬ 
creased,  but  the  error  in  the  target  dynamics  position  and 
centroid  locations  also  increased.  It  was  then  necessary  to 
increase  the  target  dynamics  velocity  variances  in  compensa¬ 
tion.  This  increase  in  the  target  dynamics  velocity  var¬ 
iance  then  resulted  in  an  increase  in  the  mismatch  between 
the  filter-computed  error  variance  and  the  actual  error 
variance.  After  several  iterations,  a  final  value  of  800 
was  selected  for  the  target  dynamics  velocity  variance  and 
0.8  for  the  jitter  variance.  The  final  tuning  run  results 
are  shown  in  Table  5.5  and  Figures  B.31-B.40.  Note  that 
there  is  a  good  match  between  filter-computed  and  actual 
error  stat'i  sties,  as  evidence  by  Figures  B.31  and  B.32.  The 
tuning  to  achieve  accuracy  in  locating  the  centroid  location 
has  also  been  preserved,  as  shown  by  the  centroid  x  and  y 
mean  errors  after  update  in  Table  5.5  and  Figures  B.39  and 
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Table  5.5.  Run  4  Temporally  Averaged  Statistics 


B.  40  . 
tuning 


Error  in  Estimate  of : 

Mean 

1  Sigma 

*(tt  ) 

• . 22712 

1 . 1871 

v(tt  ) 

.  41019 

1  .  1526 

x(t-) 

-.21384 

1 . 0945 

y(  V ) 

. 31213 

1 . 0663 

) 

. 15671 

1.0236 

> 

. 17855 

.95695 

.04210 

.21210 

.00110 

. 18875 

Comparing  the  results  in 

Table 

5.5  to  the 

run  results,  the  tuning 

process 

seems  to 

prior 

benefit 


the 


error  means  more  than  the  one  sigma  values.  At  this  point 


the  filter  was  considered  tuned  and  further  tuning  was 


discontinued . 


5 . 4  Mass  Center  Offset  Filter  Sensitivity  Study 

The  main  thrust  of  this  research  is  to  examine  the 
performance  of  a  filter  designed  to  estimate  the  offset 
between  the  missile  hardbody  center  of  mass  and  the  filter 
computed  image  centroid,  and  the  effects  on  this  performance 
of  variations  in  parameters  that  define  the  tracking  scen¬ 
ario.  Before  describing  the  results  of  these  parameter 
studies,  some  background  material  will  be  presented. 

For  this  work,  it  was  decided  to  continue  using  the  same 
size  missile  used  by  Rizzo  [14’.  The  missile  used  in  the 
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simulation  had  a  length  of  40  meters  and  a  diameter  of  3 
meters.  As  was  previously  mentioned,  it  was  decided  to 
approximate  the  shape  of  the  missile  in  the  FLIR  plane  as  a 
rectangle.  With  the  trajectory  used  in  this  work,  such  that 
the  missile  is  basically  orthogonal  to  the  line  of  sight, 
the  rectangle  dimensions  were  approximately  the  same  as  the 
missile  dimensions. 

The  next  item  of  consideration  was  the  choice  of  beam- 
width  for  the  low  power  laser  used  to  generate  the  "measur¬ 
ement"  of  the  distance  between  the  image  centroid  and  the 
missile  center  of  mass  in  the  simulation.  The  desire  to 
choose  a  beamwidth  realistic  for  the  range  used  in  the 
simulation  (approximately  2000  kilometers)  versus  the  need 
for  simplification  in  the  initial  simulation  of  this  problem 
led  to  a  tradeoff  study.  For  ease  of  implementation  of  the 
software  simulating  the  reflection  of  the  laser,  it  was 
decided  to  use  a  beamwidth  less  than  or  equal  to  the  missile 
diameter.  To  allow  beamwidths  wider  than  the  missile  dia¬ 
meter  would  have  caused  a  large  increase  in  the  amount  of 
software  necessary  for  the  beam  simulation.  It  was  also 
desired  to  have  as  large  a  number  of  simulated  reflections 
during  the  simulation  as  possible.  This  is  a  function  of 
the  accuracy  of  estimating  the  image  intensity  centroid  and 
target  velocity  (to  establish  l he  line  along  which  to  sweep 
the  laser),  the  true  offset  between  the  intensity  centroid 
and  missile  center  of  mass,  and  the  laser  beamwidth.  A 
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study  was  then  conducted  to  examine  these  tradeoffs.  Al¬ 
though  the  details  of  this  study  will  not  be  presented,  a 
value  of  2.75  meters  was  finally  selected  as  the  beamwidth 
value  at  the  target  location  to  be  used  in  the  study. 

The  sensitivity  study  examined  the  effect  of  simulta¬ 
neously  changing  both  filter-assumed  measurement  noise 
variance  and  the  truth  model  measurement  noise  variance. 

The  variance  of  the  filter  discrete-time  noise  was  set  at 
four  meters2  [8]  for  the  entire  study.  With  no  prior 
knowledge  of  the  physical  parameters  associated  with  the 
"real  world"  offset,  a  value  equal  to  ten  percent  of  the 
missile  length  was  used.  Three  measurement  error  levels 
were  examined,  representing  high,  medium,  and  low  levels  of 
accuracy.  The  measurement  error  standard  deviation  levels 
studied  were  0.2,  2,  and  20  meters.  For  this  study,  two 
meters  is  considered  the  nominal  error  [8].  The  results  of 
the  study  art  presented  in  Table  5.6  and  Figures  C.1-C.9. 

As  was  previously  mentioned  in  Chapter  4,  the  error  reported 
in  both  the  time  averaged  statistics  and  the  performance 
plots  is  the  error  between  the  truth  model  center  of  mass 
location  and  the  estimated  center  of  mass  location  which  was 
determined  using  the  filter-estimated  offset  distance,  the 
estimated  velocity  vector  orientation  angle,  and  the  es¬ 
timated  image  intensity  centroid  location.  In  Table  5.6,  the 
measurement  error  levels  are  in  meters  and  the  statistics 
are  in  pixels. 
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Table  5.6  Sensitivity  Study  Time 

Averaged 

Statistics 

Measurement  Noise 

Center  of  Mass 

Mean 

1  Sigma 

Standard  Deviation 

Error  Statistic 

0.2 

xcm  error^i  ) 

. 68248 

. 20311 

xcm  error^i  ) 

.  09571 

.18694 

2.0 

Kcm  error^i  ) 

.74503 

.22147 

xcm  error^i  ) 

. 10714 

. 14196 

20 

K cm  error  ^  i  ) 

1 . 2707 

. 39839 

xcm  error^i  ) 

.  19072 

.14516 

In  examining  the  statistics  for  the  first  two  measure¬ 
ment  precision  levels,  as  well  as  Figures  C.2,  C.3,  C.5,  and 
C.6,  several  observations  can  be  made.  First,  the  mean 
error  (particularly  after  update)  changes  only  approximately 
ten  percent  for  an  order  of  magnitude  change  in  the  measure¬ 
ment  precision  level.  Observing  Figures  C.2  and  C.4,  it  is 
apparent  that  with  increasing  measurement  noise  strength, 
the  filter  takes  longer  to  learn  the  underlying  dynamics  and 
reach  the  steady  state  mean  error  level,  as  would  be  ex¬ 
pected.  The  third  measurement  precision  level  had  much 
poorer  statistics.  In  this  case,  the  filter  is  still  learn¬ 
ing  the  underlying  dynamics  and  never  does  reach  the  steady 
state  mean  error  during  the  simulation  time  interval.  This 
is  well  illustrated  by  Figure  C.8.  Also  note  that  there  is 
little  difference  in  the  standard  deviations  between  the  two 
and  twenty  meter  measurement  precision  levels.  This  indi¬ 
cates  that  increasing  the  value  of  measurement  noise  var- 
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iance  predominantly  affects  the  amount  of  time  the  filter 
needs  to  estimate  the  offset,  and  not  the  accuracy  of  that 
estimate.  Also  note  the  effect  of  using  an  initial  filter 
covariance  of  zero,  especially  in  the  case  of  the  20  meter 
measurement  precision  level. 

5 . 5  Mass  Center  Offset  Filter  Robustness  Study 

In  a  robustness  study,  the  truth  model  is  given  para¬ 
meters  different  irom  those  in  the  filter,  and  the  perfor¬ 
mance  of  the  filter  analyzed.  This  analysis  indicates  the 
degree  of  filter  robustness  when  the  parameters  assumed  by 
the  filter  models  differ  from  the  "real"  world  as  represen¬ 
ted  by  the  truth  model.  In  this  robustness  study,  the 
measurement  error  level  was  also  the  parameter  varied. 
Filters  with  the  three  measurement  error  levels  (  measure¬ 
ment  noise  standard  deviations  of  0.2,  2,  and  20  meters) 
were  studied.  In  each  case  of  filter-assumed  measurement 
precision,  the  two  remaining  error  levels  were  consecutively 
used  in  the  truth  model,  and  the  performance  of  the  filter 
analyzed.  As  in  the  sensitivity  analysis,  the  variance  of 
the  filter  driving  noise  was  held  constant  at  four  meters 
for  all  cases  considered  in  the  robustness  study.  The 
results  are  presented  in  Table  5.7  and  Figures  D.1-D.18. 

The  parameter  values  and  statistics  from  the  sensitivity 
study  are  included  as  a  reference.  The  units  for  Table  5.7 
are  the  same  as  Table  5.6. 
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Table  5.7 


Robustness  Study  Time  Averaged  Statistics 


Measurement  Noise 
Standard  Deviation 
Filter  Truth  Model 

Mass  Center 

Error 

Statistic 

Mean  1 

Sigma 

0.2 

0.2 

xcm  error^i  ) 

. 68248 

.20311 

xcm  error ^ i  ) 

. 09571 

.18694 

0 . 2 

2.0 

cm  error'  i  1 

. 67927 

. 20276 

xcm  error (( i  ) 

.  09551 

.19553 

0 . 2 

20 

xcm  error  ( f- 1  ) 

.  78159 

.  44836 

xcm  error^i  ) 

.  10935 

.  25052 

2 . 0 

0 . 2 

xcm  error^i  ) 

. 74923 

.  22217 

x cm  error ^  i  ) 

. 10745 

.  14231 

2 . 0 

2 . 0 

x  cm  error  ( ^  ) 

.  74503 

.22147 

xcm  error^i  ) 

. 10714 

.  14196 

2.0 

20 

xcm  error i  ) 

.72303 

. 26899 

x cm  error^i  ) 

.  10539 

.22603 

20 

0 . 2 

xcm  error  ( ^  i  ) 

1 . 21810 

.40208 

x cm  error^i  ) 

.  19252 

.14554 

20 

2.0 

xcm  error^i  ) 

1 . 28010 

.  40172 

xcm  error^i  ) 

.  19235 

.  14487 

20 

20 

x cm  error  (  *  i  ) 

1 . 27070 

.  39839 

xcm  error^i  ) 

.  19072 

.  14516 

For  the  first  case,  the  trend  is  for  decreasing 
performance  with  increasing  valut.s  of  measurement  noise 
variance  in  the  truth  model.  This  is  to  be  expected  since 
the  mismatch  between  filter  and  truth  model  increases  with 
increasing  measurement  noise  variance.  However,  the  in¬ 
creases  in  the  mean  error  are  only  13  percent  for  an  order 
of  magnitude  increase  in  the  truth  model  noise  standard 
deviation.  This  is  due  to  the  fact  that  the  filter  weights 
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the  measurements  so  heavily  that,  even  with  the  increased 
uncertainty  in  the  truth  model  measurements,  the  filter  is 
quickly  able  to  learn  the  value  of  the  offset.  Note  also  in 
Figures  D.5  and  D.6  the  significant  reduction  in  the  stan¬ 
dard  deviation  after  update  for  the  case  where  20  meters  was 
used  in  the  truth  model. 

The  second  case  is  interesting  because  it  shows  the 
effect  of  having  both  larger  and  smaller  measurement  noise 
variances  in  the  truth  model  as  compared  to  the  filter. 

Here,  the  mean  errors  (both  before  and  after  update)  are 
smaller  when  the  truth  model  measurement  error  level  exceeds 
the  filter  measurement  noise  variance,  but  with  higher  error 
standard  deviations.  This  can  be  seen  in  both  in  Table  5.7 
and  Figures  D.8,  D.9,  D.ll,  and  D.12. 

The  third  case  considers  the  performance  of  a  filter 
which  has  low  confidence  in  the  measurements  it  receives. 
Note  that  the  long  transient  for  the  case  of  the  filter  with 
the  20  meter  error  standard  deviation  shows  clearly  the 
deleterious  impact  of  using  zero  as  the  initial  value  for 
the  covariance.  Here  the  mean  error  is  approximately  twice 
that  in  case  one  and  two.  Particularly  note  Figures  D.14, 
D.15,  D . 17 ,  and  D.18.  Again,  similar  to  the  third  case 
considered  in  the  sensitivity  study,  the  filter  is  still 
learning  the  value  of  the  offset  state,  and  does  not  reach 
the  steady-state  value  during  the  simulation.  However,  the 
standard  deviations  in  the  two  variations  considered  here 
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are  very  close  to  the  first  part  of  case  two.  This  situa¬ 
tion  is  similar  in  that  the  filter-assumed  measurement  noise 
variance  was  greater  than  that  in  the  truth  model. 

5 . 6  Pixel  Constant  Study 

The  next  parameter  to  be  investigated  is  the  pixel 
proportionality  constant.  The  main  purpose  of  this  study  is 
to  determine  the  FLIR  sensor  element  resolution  necessary 
for  implementation  of  a  tracking  scheme  involving  active 
illumination  of  the  target.  In  this  study,  the  variance  of 
the  filter  discrete-time  noise  was  continued  at  a  value  of 
four  meters2.  Four  meters2  (corresponding  to  the  nominal 
error  standard  deviation  of  two  meters)  was  used  as  the 
measurement  noise  variance  in  the  mass  center  offset  filter. 
The  same  value  was  also  used  for  the  measurement  noise 
variance  in  the  truth  model.  The  performance  of  the  offset 
filter  was  examined  as  the  value  of  the  pixel  proportion¬ 
ality  constant  was  varied  between  7.5,  15,  and  30  micro¬ 
radians/pixel.  The  variation  involving  the  15  micro-radian- 
/pixel  value  has  already  been  presented  in  the  sensitivity 
study,  but  is  repeated  here  to  aid  in  the  comparisons.  The 
results  are  presented  in  Table  5.8  (in  units  of  pixels)  and 
Figures  E.1-E.9. 

To  compare  the  results  in  terms  of  constant  pixel  size, 
the  statistics  in  Table  5.8  for  the  7.5  micro-radian/pixel 
case  were  determined  by  multiplying  the  mean  and  standard 
deviation  from  that  analysis  by  two.  In  a  similar  manner, 
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Table  5.8  Pixel  Constant  Study  Time  Averaged  Statistics 


Pixel  Constant 

CM  Error  Statistic 

Mean 

1  Sigma 

7.5 

A 

x  cm 

error (  ^ 

) 

2 . 57560 

1 . 27740 

xcm 

error ( * i 

) 

1 .14648 

.  34514 

15 

xcm 

A 

error ( ^ 

) 

. 74503 

.  22147 

xcm 

error ( *  i 

) 

. 10714 

.14196 

30 

*  cm 

A 

error (  ^ 

) 

. 39447 

.24665 

Kcm 

error ^  i 

) 

.01208 

.03493 

the  statistics  in  Table  5.8  for  the  30  micro-radian/pixel 
case  were  computed  by  dividing  the  mean  and  standard  devia¬ 
tion  from  that  analysis  by  two.  After  converting  to  a 
constant  pixel  si2e  (30  meters  by  30  meters,  for  the  nominal 
15  micro-radian/pixel  case  at  a  target  range  of  2000  kilo¬ 
meters),  several  observations  can  be  made.  First,  in  ob¬ 
serving  Table  5.8,  the  mean  center  of  mass  position  error 
decreased  with  increasing  pixel  size,  both  before  and  after 
update.  The  same  trend  was  true  for  the  standard  devia¬ 
tions,  except  in  comparing  the  15  and  30  micro-radian/pixel 
preupdate  standard  deviations.  These  results  seem  to  be 
contrary  to  the  expectation  of  increasing  filter  performance 
with  increasing  sensor  resolution  (i.e,  a  smaller  pixel 
constant).  However,  a  different  trend  is  apparent  in  Fi¬ 
gures  E.1-E9.  Plots  E.l  through  E.6  appear  to  echo  the 
trend  seen  in  Table  5.8  when  pixel  scale  factors  are  con¬ 
sidered.  Figures  E.7-E.9  indicate  a  much  different  trend. 
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Here  the  actual  error  is  increasing  with  time  as  shown  in 
Figure  E.7.  Figure  E.8  echoes  this,  as  the  standard  devia¬ 
tion  about  the  mean  offset  value  is  also  increasing  with 
time.  The  peaks  and  gaps  in  Figure  E.9  (caused  by  the  fact 
that  the  filter  may  go  for  extended  periods  without  receiv¬ 
ing  a  measurement  update)  would  seem  to  indicate  that,  at 
this  level  of  resolution,  the  ability  of  the  six-state 
filter  to  provide  information  about  the  location  of  the 
image  centroid  and  tne  target  velocity  orientation  angle  is 
seriously  degraded.  With  this  degradation,  the  number  of 
times  a  reflection  occurs  (a  10  Monte-Carlo  run  average  of 
8.6  reflections  for  a  seven  second  simulation,  with  210 
reflections  possible)  is  reduced,  and  the  offset  filter 
performance  suffers  accordingly. 

5 • 7  Offset  Distance  Study 

The  final  study  conducted  examined  the  mass  center 
offset  filter  performance  as  a  function  of  the  actual  offset 
distance  between  the  missile  center  of  mass  and  image  center 
of  intensity.  The  expectation  was  that,  as  offset  distance 
was  increased,  the  performance  of  the  filter  estimating  the 
offset  distance  would  suffer.  This  is  due  to  the  fact  that 
the  error  in  the  six-state  filter  estimate  of  the  velocity 
orientation  angle  used  to  aim  the  laser  increases  linearly 
with  offset  distance.  As  the  error  increases,  there  are 
fewer  reflections  which  can  be  used  to  update  the  filter. 
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The  filter  and  truth  model  noise  parameters  are 
identical  to  those  in  the  pixel  constant  study.  The  nominal 
value  of  the  pixel  constant,  15  micro-radians/pixel,  was 
also  used  in  the  study.  The  two  offset  distances  considered 
were  86  and  172  meters.  These  were  selected  as  multiples  of 
the  nominal  43  meter  value  used  for  all  the  other  analyses, 
and  because  these  values  would  allow  both  the  missile  and 
the  FLIR  image  to  be  contained  in  an  8  x  8  pixel  field  of 
view . 

The  time  averaged  statistics  are  presented  in  Table  5.9 
and  the  performance  plots  in  Figures  F.1-F.6.  In  Table  5.9, 
the  offset  distance  is  in  meters  and  the  statistics  are  in 
pixels.  Although  not  part  of  this  study,  the  statistics  for 
the  nominal  43  meter  value  are  included  in  Table  5.9  for 
reference . 

The  statistics  in  Table  5.9  confirm  the  expectation  of 
decreasing  filter  performance  with  increasing  offset 
distance.  Both  the  offset  mean  and  standard  deviation 
increase  with  increasing  offset  distance.  The  performance 
plots  also  confirm  this  idea.  In  comparing  Figures  F.1-F.3 
to  F.4-F.6,  note  that  the  filter  using  the  172  meter  offset 
is  taking  much  longer  to  reach  steady  state  performance. 

This  is  due  to  the  fact  that  there  are  fewer  measurements  (a 
ten  Monte-Carlo  run  average  of  23  reflections  per  7  second 
simulation  for  the  172  meter  offset  against  an  average  of 
28.9  for  the  86  meter  offset  case  and  32.6  for  the  43  meter 
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Table  5.9  Offset  Distance  Study  Time  Averaged  Statistics 


Offset  Distance 


Error  Statistic 


Mean 


43 


86 


172 


1  Sigma 


cm  error^i  ) 

.75403 

.  22147 

cm  error ) 

.10714 

.14196 

cm  error  (*i  ) 

.  98793 

.  24490 

cm  error  ( ^  i  ) 

.11227 

. 22204 

cm  error i  ) 

1 .61800 

.  6  6  76  6 

cm  error^i  ) 

.  12636 

. 29407 

offset  case)  for  the  filter  to  learn  the  state.  In  examin¬ 


ing  Figures  F.4  and  F.5,  note  the  large  change  in  the  fil¬ 
ter's  performance  at  approximately  three  seconds.  To  deter¬ 
mine  the  reason  for  this,  the  times  at  which  measurements 


occurred  were  examined  for  all  ten  Monte-Carlo  runs.  For 


two  of  those  runs,  only  one  measurement  was  received  by  the 
filter  in  the  first  three  seconds.  Additionally,  for  sev¬ 
eral  of  the  runs,  a  very  small  percentage  of  the  total 
number  of  measurements  in  the  simulation  occurred  before 


three  seconds.  With  only  10  runs,  these  anomalies  have  a 
significant  impact  on  the  statistics  plotted. 


5 . 8  Summary 

This  chapter  presented  the  results  of  this  thesis 
research.  The  observability  analyses  and  filter  tuning 
process  used  to  develop  the  six-state  FLIR  data  filter  were 
presented  first.  Then,  using  this  filter's  outputs  as  the 
means  to  produce  "measurements"  of  the  low-power  laser 
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returns  from  the  missile  hardbody ,  the  studies  involving  the 
mass  center  offset  filter  were  then  described.  These  stud¬ 
ies  included  a  sensitivity  and  robustness  analysis  of  the 
filter  itself,  and  studies  involving  variations  of  the  pixel 
constant  and  truth  model  offset  distance. 
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6 .  Conclusions  and  Recommendations 

6 . 1  Introduction 

This  chapter  summarizes  the  conclusions  reached  in  this 
thesis  and  recommends  topics  for  further  study.  Section  6.2 
draws  conclusions  based  on  the  analyses  of  Chapter  5. 

Section  6.3  continues  with  suggestions  for  continued 
research  in  the  area  of  FLIR  tracking  of  ballistic  missiles 
aided  by  active  low  power  laser  illumination. 

6 . 2  Conclusions 

Various  conclusions  have  been  made  in  Chapter  5.  These 
conclusions  will  now  be  brought  together  and  presented. 

6.2.1  Six-State  Filter  Tuning.  Tuning  the  filter  via 
the  target  state  dynamics  driving  noise  and  the  atmospheric 
jitter  driving  noise  proved  adequate.  The  final  tuning 
yielded  a  filter  that  could  accurately  identify  the  location 
of  the  FLIR  image  centroid.  However,  the  development  of  the 
simulation  of  the  laser  reflection  measurement  required  an 
accurate  estimate  of  both  the  image  centroid  location  and 
the  velocity  vector.  In  estimating  velocity,  the  filter's 
performance  was  less  than  adequate.  The  initial  goal  [8] 
was  to  be  able  to  aim  the  laser  accurately  enough  to  receive 
a  reflection  measurement  better  than  35  percent  of  the 
simulation  time.  The  percentage  of  measurements  actually 
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received  was  on  the  order  of  15  percent,  much  less  than 
anticipated.  Comparison  of  the  truth  model  image  centroid 
and  velocity  orientation  angle  to  the  filter  estimates  of 
these  parameters  showed  the  errors  in  the  filter-estimated 
velocity  orientation  angle  to  be  the  reason  for  the  low 
percentage  of  reflections  received.  This  problem  will  be 
discussed  further  in  Section  6.3. 

6.2.2  Offset  Filter  Sensitivity  Study.  This  study 
examined  the  effects  of  different  sensor  measurement  noise 
levels  on  filter  performance.  The  analysis  indicated  that 
the  steady-state  mean  error  in  the  offset  length  is  insensi¬ 
tive  to  measurement  noise  variance.  The  standard  deviations 
reported  were  also  relatively  independent  of  measurement 
noise  variance.  The  measurement  noise  variance  does  strong¬ 
ly  affect  the  amount  of  time  the  filter  takes  to  reach  the 
steady-state  error  level.  The  conclusion  from  this  is  that 
sensor  accuracy  may  be  a  secondary  consideration  in  system 
design  trade-off  studies,  for  the  models  used  here. 

6.2.3  Offset  Filter  Robustness  Study.  In  a  robustness 
study,  the  primary  area  of  concern  is  the  filter's  perfor¬ 
mance  when  the  "real  world",  as  represented  by  the  truth 
model,  is  different  from  the  filter's  internal  world  model. 
The  parameter  used  in  this  study  was  the  measurement  noise 
variance.  Three  filters,  with  high,  medium,  and  low  assumed 
measurement  noise  variances  were  studied.  In  each  case,  the 
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filter's  performance  was  analyzed  when  the  truth  model  was 
given  a  measurement  noise  level  different  from  that  in  the 
filter.  The  results  indicated  that  the  filter's  performance 
was  fairly  robust  to  changes  in  this  truth  model  parameter. 

These  results  lead  to  the  conclusion  that  single  filter 
performance  is  adequate  and  a  multiple  model  adaptive  filter 
algorithm  or  other  form  of  filter  that  adapts  to  this  para¬ 
meter  is  not  necessary. 

6.2.4  Pixel  Proportionality  Constant  Study.  The  re¬ 
sults  of  this  study  are  not  conclusive.  The  general  expec¬ 
tation  for  FLIR  tracking  systems  is  that  increased  sensor 
resolution  should  lead  to  increased  tracking  accuracy. 

However,  the  results  did  not  support  this  in  all  respects. 

The  intermediate  sensor  resolution  gave  the  best  filter 
performance.  No  immediate  explanation  was  found  for  this 
result.  Since  the  six-state  filter  was  tuned  using  this 
value,  it  is  unknown  whether  or  not  this  introduced  a  bias 
in  the  results.  No  conclusions  were  reached;  this  topic 
will  be  will  addressed  further  in  Section  6.3. 

6.2.5  Offset  Distance  Study.  The  final  study  in  this 
research  examined  the  effect  of  increasing  the  offset  dis¬ 
tance  between  missile  hardbody  and  infrared  intensity  cen¬ 
troid.  The  offset  distances  studied  were  selected  to  allow 
both  the  FLIR  intensity  image  and  the  reflections  from  the 
missile  hardbody  to  occur  within  an  8  x  8  pixel  field  of 
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view.  The  results  support  the  conclusion  that,  as  offset 
distance  increases,  filter  performance  decreases;  the  amount 
of  degradation  was  quantified. 

6 . 3  Recommendations 

The  recommendations  made  here  are  to  suggest  possible 
lines  of  research  to  clarify  problems  uncovered  in  this 
research  or  to  expand  preliminary  investigations  started  in 
this  work. 

6.3.1  Six-State  Filter  Modeling  and  Tuning.  Based  on 
Rizzo's  recommendations  [14],  a  reduced  order  six-state 
filter  was  finally  implemented  by  removing  acceleration 
states  frcm  the  filter  design  model.  However,  Rizzo  also 
recommended  remodeling  the  acceleration  states  using  a 
constant  turn  rate  instead  of  a  first-order  Gauss-Markov 
model.  This  suggestion  was  not  implemented  in  this  re¬ 
search.  Remodeling  the  acceleration  states,  retuning  the 
filter,  and  then  doing  an  observability  analysis,  represents 
an  alternative  approach  to  addressing  the  observability  pro¬ 
blem.  This  approach  might  also  yield  a  filter  which  better 
estimates  the  velocity  orientation  angle.  It  cannot  be 
emphasized  enough  that  accurate  velocity  vector  estimation 
is  crucial  to  the  viability  of  the  algorithm  in  it's  present 
form . 

The  second  recommendation  is  to  retune  the  six-state 
filter  to  yield  better  velocity  estimates.  The  tuning 
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accomplished  in  this  research  paid  close  attention  only  to 
the  filter's  performance  with  regard  to  estimating  the  FLIR 
intensity  image  centroid  location.  The  filter  should  be 
retuned  for  both  accurate  centroid  position  and  velocity 
vector  estimates.  For  active  illumination  of  the  target  to 
be  a  viable  concept,  both  the  image  centroid  and  the 
velocity  vector  must  be  accurately  estimated. 

Direct  examination  of  jitter  position  plots  would  also 
benefit  the  tuning  process.  In  this  effort,  the  jitter 
position  state  tuning  response  could  only  be  examined  in¬ 
directly  by  looking  at  the  image  centroid  plots.  By  direct¬ 
ly  examining  plots  of  both  dynamics  and  jitter  position 
states,  a  clearer  understanding  of  the  effect  on  filter 
tuning  of  changing  parameters  in  these  two  processes  could 
be  achieved. 

6.3.2  Illumination  Modeling  Improvement.  For  this 
thesis  research,  a  simple  geometrical  approach  was  used  to 
model  the  laser  illumination  of  the  target.  Many  effects 
were  neglected.  For  future  work,  a  more  accurate  model  is 
needed.  This  model  should  include  more  accurate 
representation  of  the  missile's  three-dimensional  shape. 

The  simulation  of  the  laser  reflection  itself  should  include 
missile  shape  effects  (laser  reflecting  off  cylindrical 
object  rather  than  flat  surface),  and  the  effects  of  atmo¬ 
spheric  distortion  in  the  reflected  laser  light.  Also,  the 
dynamics  of  sweeping  the  laser  up  and  down  the  estimated 
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velocity  vector  direction  from  the  intensity  centroid  should 
be  incorporated.  Alternate  search  patterns,  in  which  the 
laser  is  swept  both  along  and  perpendicular  to  the  estimated 
velocity  vector  in  a  sinusoidal  pattern,  should  be  con¬ 
sidered.  This  would  increase  the  probability  of  reflection 
by  increasing  the  search  area,  as  well  as  reducing  the 
dependence  of  the  algorithm  on  an  accurate  velocity  vector 
estimate.  The  increased  accuracy  of  illumination  modeling 
will  significantly  enhance  the  credibility  of  the  simulation 
and  conclusions  drawn  from  it. 

6.3.3  Different  Initial  Filter  Parameters.  With  no 
previous  knowledge  of  the  initial  errors  involved  with 
pointing  a  laser  at  a  missile,  the  initial  value  for  the 
offset  filter  mean  and  covariance  were  set  to  zero.  As  was 
evident  in  the  sensitivity  and  robustness  studies,  this  had 
a  significant  effect  on  filter  performance.  Based  on  the 
steady  state  filter  errors  found  in  this  research,  it  is 
recommended  that  the  analyses  be  performed  again  using  an 
initial  mean  of  0.1  pixels  and  covariance  of  .05  pixels'  to 
confirm  the  results  found  here. 

6.3.4  Pixel  Proportionality  Constant  Filter  Tuning. 

The  results  of  the  pixel  constant  study  were  inconclusive. 

To  resolve  the  issue  of  sensor  resolution,  the  six-state 
filter  should  be  separately  tuned  using  the  two  pixel  con¬ 
stants  studied  in  Section  5.6.  After  tuning,  comparison  of 
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the  performance  of  the  three  filters  should  help  resolve 
this  issue. 

6.3.5  Mass  Center  Error  Resolution.  The  error  between 
the  filter  estimated  center  of  mass  and  the  truth  model 
center  of  mass  used  in  the  analysis  was  calculated  as  the 
scalar  length  between  two  coplanar  points.  To  aid  in  future 
work,  this  error  should  be  resolved  into  a  component  along 
the  estimated  velocity  vector,  and  a  component  perpendicular 
to  the  velocity  vector.  By  doing  this,  it  should  be  pos¬ 
sible  to  determine  the  contribution  of  the  offset  distance 
estimate  (along  velocity  vector  component)  and  the  estimated 
velocity  vector  orientation  angle  (perpendicular  to  velocity 
vector  component)  to  the  total  error. 
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Appendix  A.  Mass  Center  Measurement  Simulation 


As  presented  in  Section  2.4.2,  the  mass  center  measure¬ 
ment  simulation  is  developed  in  terms  of  plane  geometry. 

The  projection  of  the  missile  hardbody  shape  onto  the  two- 
dimensional  FLIR  plane  is  represented  by  a  rectangle.  To 
simulate  the  reflection  of  the  laser,  and  the  mass  center 
measurement  derived  from  that  reflection,  the  path  the  laser 
takes  as  it  is  dithered  along  the  velocity  vector  had  to  be 
simulated.  It  was  decided  to  represent  the  laser's  path  as 
a  rectangle  in  the  FLIR  plane  also.  The  width  of  the  rec¬ 
tangle  was  chosen  to  approximate  the  laser  beamwidth  at  the 
target.  For  this  study,  a  beamwidth  of  2.75  meters  was 
used.  The  length  of  the  rectangle  was  determined  using 
three  times  the  truth  model  offset  distance  between  the 
missile  center  of  mass  and  the  FLIR  image  centroid.  This 
was  done  to  insure  that  the  rectangle  was  long  enough  for 
intersection  with  the  rectangle  representing  the  missile  to 
occur . 

The  rectangle  representing  the  laser  path  is  located  on 
the  FLIR  plane  using  the  filter  estimated  image  centroid 
location  and  velocity  vector  orientation.  This  is  shown  in 
Figure  A.l.  One  end  of  the  rectangle  is  located  coincident 
with  the  image  centroid  location.  The  coordinates  of  the 
opposite  end  of  the  laser  path  rectangle  are  located  using 
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the  filter  estimated  velocity  vector  orientation  angle.  The 
calculations  for  the  endpoint  are  as  follows: 

*  endpoint  =  * centroid  +  len*th  cos  9/ 
y endpoint  ~  ^ centroid  ^enS^  sin  8^ 

where : 

K  endpoint '  V endpoint  =  the  FLIR  Plane  coordinates  of  the 

rectangle  endpoints 

* centroid '  V centroid  =  the  FLIR  Plane  ima^e  centroid 

coordinates 

length  =  desired  rectangle  length,  in  pixels 
=  filter  estimated  velocity  vector 
orientation  angle 

Once  the  endpoints  have  been  determined,  the  coordinates  of 
the  corners  of  the  rectangle  are  determined  by  offsetting 
half  the  beamwidth  distance  from  the  endpoints  along  a  line 
perpendicular  to  the  velocity  vector. 

The  simulation  of  the  laser  reflection  and  the  correspo¬ 
nding  center  of  mass  measurement  is  done  in  several  steps. 
First,  the  laser  beam  striking  the  missile  and  reflecting  is 
simulated  by  determining  the  area  of  intersection  between 
the  missile  rectangle  and  the  rectangle  representing  the 
laser  beam  path,  as  shown  in  Figure  A. 2.  This  area  of 
intersection  is  determined  using  plane  geometry.  Figure  A. 3 


(A-l) 
( A- 2  ) 
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Figure  A. 3.  Center  of  Mass  Measurement  Simulation 
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shows  an  expanded  view  of  the  area  of  intersection.  Note 
that,  for  this  example,  two  of  the  corners  of  the  perimeter 
of  the  intersection  area  are  determined  by  the  intersections 
of  line  segments  representing  the  outline  of  the  missile  and 
laser  path  rectangle.  The  determination  of  the  intersection 
coordinates  is  a  matter  of  plane  geometry.  For  example,  the 
determination  of  the  x  and  y  coordinates  of  point  1  is  made 
as  follows: 

The  slope  of  the  line  segment  between  (x.  ,y.  )  and  point  two 
is  m,  ,  where  m,  =  tan  6^;  the  equation  of  this  line  is 

y  =  m  (x-xL)  +  y.  ( A-3  ) 

Likewise,  the  slope  of  the  line  segment  between  ( x_ , y  )  and 


point  three  is  m2 ,  where  m£ 
this  line  is 


_ 1 _ . 

tan  ' 


and  the  equation  of 


y  =  m2(x-xj  +  y_  (A-4) 

Setting  Equation  (A-3)  equal  to  Equation  (A-4)  and  solving 
for  x  yields: 

m  x.  -  mx_  +  y  -  y 

, ,  T*  -  /  X  _  (Z  \ 


Solving  Equation  (A-5)  and  substituting  this  value  into 
either  Equation  (A-3)  or  (A-4)  yields  the  coordinates  of 
point  one  in  Figure  A. 3.  The  coordinates  for  point  two,  or 
for  any  other  intersection  point,  are  determined  in  a  simi¬ 
lar  manner.  One  area  of  concern  is  that  Equations  (A-3)-(A- 
5)  involve  the  equations  of  lines,  not  line  segments. 
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Therefore,  in  the  simulation,  the  software  checks  the  value 
of  the  x  intercept  point  to  determine  if  it  is  between  the  x 
coordinate  of  point  three  and  x*.  If  this  condition  is 
true,  then  the  line  segments  do  intersect,  some  portion  of 
the  laser  beam  rectangle  falls  on  the  missile  rectangle,  and 
reflection  occurs .  If  this  condition  is  not  met,  then  the 
intersection  point  occurs  beyond  the  line  segment  endpoints. 
In  this  case,  no  portion  of  the  laser  rectangle  and  the 
missile  rectangle  are  in  common.  This  constitutes  a  "miss", 
and  no  reflection  occurs. 

Assuming  reflection  has  occurred,  the  next  step  is  to 
determine  the  center  of  mass.  Since  the  entire  purpose  of 
using  low  power  laser  reflections  is  to  obtain  information 
about  the  center  of  mass  of  the  missile,  the  center  of  the 
intersection  area  is  eguated  to  the  "measured"  center  of 
mass.  For  the  example  of  Figure  A. 3,  the  area  defined  by 
points  one,  two,  and  three  is  triangular.  Approximating  the 
center  of  this  area  is  basically  approximating  the  center  of 
a  triangle.  The  center  of  a  triangle  would  normally  be 
determined  by  the  intersection  of  three  lines,  each  line 
running  from  a  corner  to  the  opposite  side,  and  that  line 
bisecting  the  angle  at  that  corner.  This  intersection  is 
approximated  by  first  taking  a  point  one  third  the  distance 
from  point  one  to  two  along  that  line  segment.  At  that 
point,  a  perpendicular  line  segment  is  run  to  the  line 
segment  between  points  two  and  three.  The  center  is 
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approximated  by  a  point  halfway  along  the  perpendicular  line 
segment.  The  distance  between  this  center  point  and  the 
FLIR  image  centroid  of  Figure  A.l  becomes  the  offset  length 
measurement  of  Equation  (2-44).  If  the  area  of  intersection 
is  rectangular  or  trapezoidal,  then  the  calculations  for  the 
center  of  the  are  made  accordingly. 
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Appendix  B.  Plots  for  Six-State  Filter  Tunin 
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Figure  B.7  Tuning  Run  #1  X-Minus  Centroid  Position  Error 
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Figure  B.12  Tuning  Run  #2  Y-Channel  Filter  Vs  Actual  Error 
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Appendix  C.  Plots  for  Sensitivity  Stud 
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Figure  C.6  2  Meter  Error  Offset-Plus  Error 
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Figure  C.9  20  Meter  Error  Offset-Plus  Error 


Appendix  D.  Plots  for  Robustness  Study 
Discussion  in  Section  5.5 
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Figures  D.10-D.12 

Filter  Error  20  meters,  Truth  Model  Error  .2  meters 
Figures  D.13-D.15 

Filter  Error  20  meters.  Truth  Model  Error  2  meters 
Figures  D.16-D.18 
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Figure  D.l  Offset  Filter  Vs  Actual  Error 
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Figure  D.3 


Offset  Plus  Error 
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Figure  D.9  Offset  Plus  Error 
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Figure  D.13  Offset  Filter  Vs  Actual  Error 
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Figure  D.14  Offset  Minus  Error 
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Figure  D.15  Offset  Plus  Error 
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Figure  D.16  Offset  Filter  Vs  Actual  Error 
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Figure  D.17  Offset  Minus  Error 
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Appendix  E.  Plots  for  Pixel  Constant  Stud 


Discussion  in  Section  0.6 
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Appendix  F.  Plots  for  Offset  Distance  Study 
Discussion  in  Section  5.7 

Offset  Distance  (Meters) 
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Figures  F.1-F.3 

172 

Figures  F . 4-F . 6 
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